​x

<font size="+3"><strong>Time Series: Statistical Models</strong></font>

# Autoregression

warnings.simplefilter(action="ignore", category=FutureWarning)

Since we'll be working with the `"air-quality"` data again, we need to connect to the server, start our client, and grab the data we need.

client = MongoClient(host="localhost", port=27017)

<font size="+1">Practice</font>

lagos_prac = ...

In order to get our data into a form we can use to build our model, we're going to need to transform it in several key ways. The first thing we need to do is to get the data we need, and save the results in a DataFrame. Since we're interested in predicting the changes in air quality over time, let's set the DataFrame's index to `"timestamp"`:

    {"metadata.site": 3, "metadata.measurement": "P2"},

<font size="+1">Practice</font>

​x

## Localizing the Timezone

df.index = df.index.tz_localize("UTC").tz_convert("Africa/Lagos")

## Resampling Data

df = df["P2"].resample("1H").mean().fillna(method="ffill").to_frame()

Notice the second half of the code:

# In `shift(1), the `1` is the lagged interval.

Finally, let's drop our null values:

y = df["P2"].resample("1H").mean().fillna(method="ffill")

<font size="+1">Practice</font>

df_prac["P2.L1"] = ...

## Exploring the Data

### Time Series Line Plot

Example of

### Creating an ACF Plot

fig1, ax = plt.subplots(figsize=(15, 6))

Each of the dots on our plot represents a correlation coefficient. The first data point in the top left of the graph tells us that at time-step 0, the correlation coefficient was 1, meaning that there was a perfect correlation. That makes sense, because you can't lag from time-step 0, so the coefficient can't be anything other than 1. But, starting at hour 1, the coefficient drops precipitously, and we see our autocorrelation coefficients slowly decay over time. As our lag recedes further into the past, the correlations break down; a prediction you made five hours ago about what's happening right now is going to be a lot more reliable than a prediction you made 96 hours ago.

fig2, ax = ...

### Creating a PACF Plot

fig1, ax = plt.subplots(figsize=(15, 6))

Aha! This looks very different. There are two things to notice here:

fig2, ax = ...

## Working with Rolling Windows

Rolling window statistics are very helpful in smoothing noisy data when making time series predictions. Let's see it with an example. Since we're interested in making predictions about the air quality in Lagos, it would be helpful to understand the rolling average for the PM 2.5 readings with a line plot. To keep things manageable, we'll set our window-size to one week.

# `168` is the number of hours in a week.

Even though there are lots of peaks and valleys here, we're starting to see an emerging trend.

​x

Besides rolling sum and rolling average, rolling window statistics can be applied to a lot of other statistics depends on the problem you are facing. In the example below, when we use **GARCH** model to analyze stock prices, we can use rolling window to calculate standard deviation.

### Splitting the Data in pandas

cutoff_test1 = int(len(y) * 0.95)

<font size="+1">Practice</font>

cutoff_test2 = ...

## Building the Model

### Baseline

mae_baseline = mean_absolute_error(y_train, y_pred_baseline)

<font size="+1">Practice</font>

y_prac_pred_baseline = ...

### Iterating

model = AutoReg(y_train, lags=24, old_names=False).fit()

Notice that, unlike our linear regression model which we built using scikit-learn, we're combining instantiation and fitting into one step; statsmodels includes that ability in its `AutoReg` method.

model_prac = ...

Autoregression models need us to generate **in-sample predictions** in order to calculate the MAE of our training data. In-sample predictions use data that's already part of our sample. That's to distinguish it from out-of-sample predictions, which we'll talk about a bit later. The statsmodels library includes a method called [`predict`](https://www.statsmodels.org/stable/examples/notebooks/generated/predict.html) that can help us here. Above, the `AutoReg` method includes this line:

y_pred = model.predict().dropna()

Once we've done that, we can calculate the MAE of the predictions in our training set.

training_mae = mean_absolute_error(y_train.loc[y_pred.index], y_pred)

<font size="+1">Practice</font>

    y_prac_train.loc[y_prac_pred.index], y_prac_pred

### Residuals

y_train_resid = y_train - y_pred

Now we can make a line plot of our model's residuals.

fig1, ax = plt.subplots(figsize=(15, 6))

The ideal residual plot has a random set of datapoints spread evenly on both sides of the line. The plot we just made actually looks pretty good; there are some significant outliers, but, on the whole, the bars describe an even band of values, which is what we're looking for.

<font size="+1">Practice</font>

y_prac_train_resid = ...

Let's also take a look at a histogram of the residuals to help us see how they're distributed.

y_train_resid.hist();

Remember, when we make histograms, we're trying to answer two questions:

### ACF Plots

plot_acf(y_train_resid.dropna(), ax=ax);

At first, this might seem wrong, but we're actually looking for a mostly-flat graph here. This is an indication that our model describes all the **seasonality**, or regular changes, in our data. In other words, this graph is exactly what we're looking for.

​x

​x

## Evaluating the Model

### Out-of-Sample Predictions

y_pred_test = model.predict(y_test.index.min(), y_test.index.max())

Now that we have a prediction, we can calculate the MAE of our out-of-sample data.

test_mae = mean_absolute_error(y_test, y_pred_test)

<font size="+1">Practice</font>

    y_prac_test.index.min(), y_prac_test.index.max()

Now that we have some out-of-sample predictions, we can compare it to our in-sample predictions using a line plot. The first step there is to create a new DataFrame called `test1_predictions` with two columns: one for the `y_test` data (the true data) and one for the `y_pred` (the predicted data). It's always a good idea to print the first five rows of a new DataFrame, just to make sure it looks the way it should.

    {"y_test": y_test, "y_pred": y_pred_test}, index=y_test.index

That looks correct, so we can move on to our line plot.

fig = px.line(test1_predictions, labels={"value": "P2"})

This looks kind of strange, but it's actually exactly what we would expect to see. At the beginning, the `y_pred` data has a fair amount of predictive power, but, as time goes on, the predictions become less and less accurate. It's kind of like what happened with our ACF plots, only in reverse. Last time, the model lost its predictive power as the lag increased. Here, the model loses its predictive power as the horizon — how far away from the present your predictions are — increases.  But don't worry! We'll fix it in a second.

​x

​x

### Walk-forward Validation

# Then, we define a variable that takes into account what's happened in the past

You'll notice that we're using the same `AutoReg` method we used before, only this time, we're using the `y_train` data. Also like before, the `24` is telling the model how many hours it should pay attention to. If you change that number, the MAE will change too.

print("Test MAE 1 (walk forward validation):", round(test1_mae, 2))

<font size="+1">Practice</font>

y_prac_pred_wfv = ...

test2_mae = ...

## Communicating the Results

print(model.params)

That looks pretty good, but showing it in a line plot would be much better.

    {"y_test": y_test, "y_pred": y_pred_wfv}, index=y_test.index

That looks much better! Now our predictions are actually tracking the `test` data, just like they did in the linear regression model.

​x

​x

# ARMA Models & Hyperparameters

**ARMA** stands for Auto Regressive Moving Average, and it's a special kind of **time-series** analysis. So far, we've used autoregression (AR) to build our time-series models, and you might recall that AR models rely on values that remain relatively stable over time. That is, they can predict the future very well, as long as the future looks roughly the same as the past. The trouble with predicting the future is that things can suddenly change, and as a result, the future doesn't look much like the past anymore. These sudden changes — economists call them *endogenous shocks* — can be as big as a hurricane destroying a city or an unexpected increase in the minimum wage, and they can be as small as a new restaurant opening in a neighborhood or a single person losing their job. In our data, the air quality might be changed if there was a nearby forest fire, or if a building collapsed near one of the sensors and raised a giant cloud of dust.

## Cleaning the Data

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

And then we need to get our database client up and running.

client = MongoClient(host="localhost", port=27017)

Then, we need to clean our data. All the examples will use data from Site 3; all the practice sets will use Site 2. If you need a refresher on how all those methods work, refer back to the Autoregression notebook.

df.index = df.index.tz_localize("UTC").tz_convert("Africa/Lagos")

<font size="+1">Practice</font>

df_prac["P2.L1"] = ...

## Exploring the Data

y.hist();

<font size="+1">Practice</font>

​x

This is what the data looks like when our sample is 1-hour intervals, but we might want to be able to quickly change our sample to other intervals of time. First, we'll create a function called `wrangle`, and then add an **argument**. In Python, arguments tell the function what to do. This function already has an argument called `collection`, so we'll need to add another to make resampling work. We'll call that argument `resamp_pd`. <span style='color: transparent; font-size:1%'>WQU WorldQuant University Applied Data Science Lab QQQQ</span>

# Here's where the new argument goes. We're setting the default value to `"1H"`.

Now let's change `"1H"` to `"1D"` and see what happens.

y = wrangle(lagos, resamp_pd="1D")

As you can see on the left side of the table, the samples are now at one day intervals, which is exactly what we wanted!

y.hist();

This looks pretty different! It's always nice to have a diversified dataset.

print(y_prac)

​x

Like with our AR model, we need to create ACF and PACF plots to see what's happening with the correlation coefficients.

fig, ax = plt.subplots(figsize=(15, 6))

And now let's make a PACF plot.

fig, ax = plt.subplots(figsize=(15, 6))

<font size="+1">Practice</font>

​x

​x

## Splitting the Data

# Notice that the date format is `YYYY-MM-DD`

<font size="+1">Practice</font>

y_prac_train = ...

## Building the Model

### Baseline

mae_baseline = mean_absolute_error(y_train, y_pred_baseline)

<font size="+1">Practice</font>

y_prac_pred_baseline = ...

### Iterating

### Hyperparameters

p_params = range(0, 25, 8)

And let's set our `q` values to include values from 0 to 3, moving in steps of 1:

q_params = range(0, 3)

<font size="+1">Practice</font>

p_params_prac = ...

In order to tell the model to keep going through all the possible combinations, we'll add in a pair of `for` loops. (If you need a refresher on `for` loops, refer to Notebook 001.)

        # Here's where we tell the model how we want it to deal with time

<font size="+1">Practice</font>

​x

Now that we have a working ARMA model, let's turn the output into a DataFrame so we can see what happened.

mae_grid = pd.DataFrame(maes)

And let's visualize the DataFrame in a heatmap. (If you need a refresher on how to create a heatmap in seaborn, refer to Notebook 008.)

plt.title("Grid Search (Criterion: MAE)");

<font size="+1">Practice</font>

mae_grid_prac = ...

It looks like our MAE values are in the right place, but let's try some other ways to explore our new model using the `plot_diagnostics` method.

fig, ax = plt.subplots(figsize=(15, 12))

As usual, we have quite a lot to sift through here. The first graph is showing us our model's residuals. Ideally, we'd like to see this be as close to zero as possible, and this graph is telling us that, for the most part, we have a good model.

​x

## Communicating the Results

    {"y_train": y_train, "y_pred": y_train_pred}, index=y_train.index

<font size="+1">Practice</font>

​x

# ARCH and GARCH Models

**ARCH** stands for autoregressive conditionally heteroscedastic, which models the variance of a time series. **ARCH** model assumes variance at time $t$ depends on the **past squared observations**. A **GARCH** (generalized autoregressive conditionally heteroscedastic) model uses values of the **past squared observations** and **past variances** to model the variance at time $t$. **ARCH** and **GARCH** models are widely used in finance and econometric time series analysis. For more details, check out these YouTube videos.

YouTubeVideo("Li95a2biFCU")

YouTubeVideo("inoBpq1UEn4")

Let's see an example of using **GARCH** model in forecasting Apple stock price volatility. We first import the data.

df = pd.read_csv("./data/AAPL.csv", parse_dates=["date"], index_col="date")

The dataset range from 1999 to 2022. For simplicity, we only take 5 years of data.

df = df["2015-10-13":]

## Calculating Returns

The next step is to calculate the volatility of the stock close prices. Here we can use the `pct_change` function from pandas to calculate the daily percentage change, then multiply the result by 100 to get the returns:

df["return"] = df["close"].pct_change() * 100

We can then take out the return column as our training data. Note the first observation is missing since there is no past value to observe. We need to drop it:

AAPL_return = df["return"].dropna()

We can check the histogram to see the distribution of the returns over the past five years:

plt.title("Distribution of AAPL Daily Returns");

There's a negative outlier in this date range, with the `idxmin` function, we find out it was in 31 August 2020. The stock price has a huge drop due to a [stock split](https://www.cnbc.com/2020/08/31/history-of-apple-stock-splits-says-dont-rush-in-to-buy-cheaper-shares.html).

AAPL_return.idxmin(), AAPL_return.min()

We can also check the standard deviation of the whole dataset with the pandas `std()` function:

AAPL_return.std()

To see the statistic more clearly, we can use some time series plots. In addition to the Apple stock price return plot, we can also plot rolling volatility of the return to smooth out the noise, and see how return and volatility are associated with each other. First we can calculate a 50-day rolling standard deviation series for Apple stock price return:

AAPL_return_rolling_50d_volatility = AAPL_return.rolling(window=50).std().dropna()

We can plot these two series:

    ax=ax, label="50d rolling volatility", linewidth=3

Here we can see that volatility goes up when the returns change drastically up or down. For instance, when return dropped enormously in 2020 August, volatility also increased a lot. However, this plot reveals a problem. We want to use returns to see if high volatility on one day is associated with high volatility on the following day. But high volatility is caused by large changes in returns, which can be either positive or negative. How can we assess negative and positive numbers together without them canceling each other out? The common solution used in building ARCH or GARCH models are to square the returns. Let's plot the squared returns:

(AAPL_return**2).plot(ax=ax, label="daily return")

 Now we it much easier to see groups high and low volatility and build a GARCH model. To build a GARCH model, we need to figure out both the `p` and `q` parameters. The `p` parameter is handling correlations at prior time steps and the `q` parameter deals with prior variances, like shocks. It also uses the notion of lag. To see how many lags we should have in our model, we should create an ACF and PACF plot — but using the squared returns.

fig, ax = plt.subplots(figsize=(15, 6))

fig, ax = plt.subplots(figsize=(15, 6))

Both the ACF and PACF graph show one lag is enough to build the model.

## Building the Model

model = arch_model(AAPL_return, p=1, q=1, rescale=False).fit(disp=0)

## Common Metrics

The model summary provides a lot of information about the model, like the trained parameters, model performance, etc. **AIC** and **BIC** are important measurements for model performance. **AIC** stands for Akaike Information Criterion, it measures the goodness of fit of any estimated statistical model. **BIC** is Bayesian Information Criteria that selects model from a finite set of models. When fitting models, it is possible to increase model performance by adding more parameters, which would also result in overfitting. The BIC resolves this problem by introducing a penalty term for the number of parameters in the model. The penalty term is larger in BIC than in AIC. Though BIC is always higher than AIC, lower the value of these two measures, better the model.

## Standardized Residuals

After fitting the model with the data, we can get also get the **residuals** to see how the model performs. The residual at time $t$ is the observed return at time $t$ minus the model's estimated mean return at time $t$. For other models that model observations directly, like stock prices, the assumptions are the residuals are the noises that follow a normal distribution. In GARCH model, since we are model returns, which is the variance of the stock prices, the residuals of the variance will not follow a normal distribution. So we need to use **standardized residuals** instead of residuals to check the normality. Standardized residual at time $t$ is the residual at time $t$ divided by the square root of the model estimated variance at time $t$. Let's check the plot of the standardized residuals across the time series:

model.std_resid.plot(ax=ax, label="Standardized Residuals")

The standardized residuals are moving around 0 except for the outlier events from 2020 August. Let's check the normality by the histogram:

plt.title("Distribution of Standardized Residuals");

If we exclude the outlier, the other standardized residuals do follow a normal distribution with mean 0.

## Evaluation

(-2 * model.conditional_volatility.rename("")).plot(ax=ax, color="C1")

We can then select the test set of the data and do a walk-forward validation on the GARCH model:

    next_pred = model.forecast(horizon=1, reindex=False).variance.values[0][0] ** 0.5

We can plot the predicted volatility versus return in this test set, and see the model has captured most of the variance here:

(2 * y_test_wfv).plot(ax=ax, c="C1", label="2 SD Predicted Volatility")

## Forecasting

The last step is to make the forecasts with our trained model. Let's make a one-day forecast first to see how GARCH model forecasting works:

one_day_forecast = model.forecast(horizon=1, reindex=False).variance

There are two things we need to keep in mind here. First, our model forecast shows the predicted variance, not the standard deviation / volatility. So we'll need to take the square root of the value. Second, the prediction's index is in the format of `"h.1"` as the next date of the last training data date. It's not in the form of DatetimeIndex that we desire. So we need to format the forecasts to be more readable.

We first generate 5 days predictions with our trained GARCH model:

prediction_dates = pd.bdate_range(start=start, periods=prediction.shape[1])

Then we define a function to clean the predictions to be more readable:

prediction_formatted = pd.Series(data, index=prediction_index)

# References & Further Reading

---

​x

<font size="+3"><strong>3.1. Wrangling Data with MongoDB 🇰🇪</strong></font>

from IPython.display import VimeoVideo

VimeoVideo("665412094", h="8334dfab2e", width=600)

VimeoVideo("665412135", h="dcff7ab83a", width=600)

**Task 3.1.1:** Instantiate a `PrettyPrinter`, and assign it to the variable `pp`.

pp = PrettyPrinter(indent = 2)

# Prepare Data

## Connect

VimeoVideo("665412155", h="1ca0dd03d0", width=600)

**Task 3.1.2:** Create a client that connects to the database running at `localhost` on port `27017`.

client = MongoClient(host="localhost",port=27017)

## Explore

VimeoVideo("665412176", h="6fea7c6346", width=600)

**Task 3.1.3:** Print a list of the databases available on `client`.

pp.pprint(list(client.list_databases()))

[ {'empty': False, 'name': 'admin', 'sizeOnDisk': 40960},
  {'empty': False, 'name': 'air-quality', 'sizeOnDisk': 7000064},
  {'empty': False, 'name': 'config', 'sizeOnDisk': 110592},
  {'empty': False, 'name': 'local', 'sizeOnDisk': 73728},
  {'empty': False, 'name': 'wqu-abtest', 'sizeOnDisk': 585728}]

VimeoVideo("665412216", h="7d4027dc33", width=600)

**Task 3.1.4:** Assign the `"air-quality"` database to the variable `db`.

db = client["air-quality"]

VimeoVideo("665412231", h="89c546b00f", width=600)

**Task 3.1.5:** Use the [`list_collections`](https://pymongo.readthedocs.io/en/stable/api/pymongo/database.html?highlight=list_collections#pymongo.database.Database.list_collections) method to print a list of the collections available in `db`.

for c in db.list_collections():

system.views
lagos
system.buckets.lagos
nairobi
system.buckets.nairobi
dar-es-salaam
system.buckets.dar-es-salaam

VimeoVideo("665412252", h="bff2abbdc0", width=600)

**Task 3.1.6:** Assign the `"nairobi"` collection in `db` to the variable name `nairobi`.

nairobi = db["nairobi"]

VimeoVideo("665412270", h="e4a5f5c84b", width=600)

**Task 3.1.7:** Use the [`count_documents`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.count_documents) method to see how many documents are in the `nairobi` collection.

nairobi.count_documents({})

202212

VimeoVideo("665412279", h="c2315f3be1", width=600)

**Task 3.1.8:** Use the [`find_one`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.find_one) method to retrieve one document from the `nairobi` collection, and assign it to the variable name `result`.

result = nairobi.find_one({})

{ 'P1': 39.67,
  '_id': ObjectId('6334b0e98c51459f9b198d27'),
  'metadata': { 'lat': -1.3,
                'lon': 36.785,
                'measurement': 'P1',
                'sensor_id': 57,
                'sensor_type': 'SDS011',
                'site': 29},
  'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}

VimeoVideo("665412306", h="e1e913dfd1", width=600)

**Task 3.1.9:** Use the [`distinct`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.distinct) method to determine how many sensor sites are included in the `nairobi` collection.

nairobi.distinct("metadata.site")

[6, 29]

VimeoVideo("665412322", h="4776c6d548", width=600)

**Task 3.1.10:** Use the [`count_documents`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.count_documents) method to determine how many readings there are for each site in the `nairobi` collection.

print("Documents from site 29:",nairobi.count_documents({"metadata.site":29}))

Documents from site 6: 70360
Documents from site 29: 131852

VimeoVideo("665412344", h="d2354584cd", width=600)

**Task 3.1.11:** Use the [`aggregate`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.aggregate) method to determine how many readings there are for each site in the `nairobi` collection.

        {"$group":{"_id":"$metadata.site","count":{"$count":{}}}}

[{'_id': 6, 'count': 70360}, {'_id': 29, 'count': 131852}]

VimeoVideo("665412372", h="565122c9cc", width=600)

**Task 3.1.12:** Use the [`distinct`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.distinct) method to determine how many types of measurements have been taken in the `nairobi` collection.

nairobi.distinct("metadata.measurement")

['P2', 'humidity', 'P1', 'temperature']

VimeoVideo("665412380", h="f7f7a39bb3", width=600)

**Task 3.1.13:** Use the [`find`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.find) method to retrieve the PM 2.5 readings from all sites. Be sure to limit your results to 3 records only.

result = nairobi.find({"metadata.measurement":"P2"}).limit(3)

[ { 'P2': 34.43,
    '_id': ObjectId('6334b0ea8c51459f9b1a0db2'),
    'metadata': { 'lat': -1.3,
                  'lon': 36.785,
                  'measurement': 'P2',
                  'sensor_id': 57,
                  'sensor_type': 'SDS011',
                  'site': 29},
    'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)},
  { 'P2': 30.53,
    '_id': ObjectId('6334b0ea8c51459f9b1a0db3'),
    'metadata': { 'lat': -1.3,
                  'lon': 36.785,
                  'measurement': 'P2',
                  'sensor_id': 57,
                  'sensor_type': 'SDS011',
                  'site': 29},
    'timestamp': datetime.datetime(2018, 9, 1, 0, 5, 3, 941000)},
  { 'P2': 22.8,
    '_id': ObjectId('6334b0ea8c51459f9b1a0db4'),
    'metadata': { 'lat': -1.3,
                  'lon': 36.785,
                  'measurement': 'P2',
                  'sensor_id': 57,
                  'sensor_type': 'SDS011',
                  'site': 29},
    'timestamp': datetime.datetime(2018, 9, 1, 0, 10, 4, 374000)}]

VimeoVideo("665412389", h="8976ea3090", width=600)

**Task 3.1.14:** Use the [`aggregate`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.aggregate) method to calculate how many readings there are for each type (`"humidity"`, `"temperature"`, `"P2"`, and `"P1"`) in site `6`.

        {"$group":{"_id":"$metadata.measurement","count":{"$count":{}}}}

[ {'_id': 'P2', 'count': 18169},
  {'_id': 'humidity', 'count': 17011},
  {'_id': 'P1', 'count': 18169},
  {'_id': 'temperature', 'count': 17011}]

VimeoVideo("665412418", h="0c4b125254", width=600)

**Task 3.1.15:** Use the [`aggregate`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.aggregate) method to calculate how many readings there are for each type (`"humidity"`, `"temperature"`, `"P2"`, and `"P1"`) in site `29`.

        {"$group":{"_id":"$metadata.measurement","count":{"$count":{}}}}

[ {'_id': 'P2', 'count': 32907},
  {'_id': 'humidity', 'count': 33019},
  {'_id': 'P1', 'count': 32907},
  {'_id': 'temperature', 'count': 33019}]
[]

## Import

VimeoVideo("665412437", h="7a436c7e7e", width=600)

**Task 3.1.16:** Use the [`find`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.find) method to retrieve the PM 2.5 readings from site `29`. Be sure to limit your results to 3 records only. Since we won't need the metadata for our model, use the `projection` argument to limit the results to the `"P2"` and `"timestamp"` keys only.

    {"metadata.site":29,"metadata.measurement":"P2"},

VimeoVideo("665412442", h="494636d1ea", width=600)

**Task 3.1.17:** Read records from your `result` into the DataFrame `df`. Be sure to set the index to `"timestamp"`.

df = pd.DataFrame(result).set_index("timestamp")

assert isinstance(df.index, pd.DatetimeIndex), "`df` should have a `DatetimeIndex`."

---

​x

<font size="+3"><strong>Databases: PyMongo</strong></font>

# Working with PyMongo

from pprint import PrettyPrinter

## Databases

## Servers and Clients

client = MongoClient(host="localhost", port=27017)

## Semi-structured Data

## Exploring a Database

pp.pprint(list(client.list_databases()))

[ {'empty': False, 'name': 'admin', 'sizeOnDisk': 40960},
  {'empty': False, 'name': 'air-quality', 'sizeOnDisk': 7000064},
  {'empty': False, 'name': 'config', 'sizeOnDisk': 98304},
  {'empty': False, 'name': 'local', 'sizeOnDisk': 73728},
  {'empty': False, 'name': 'wqu-abtest', 'sizeOnDisk': 585728}]

It looks like this server contains four databases: `"admin"`, `"air-quality"`, `"config"`, and `"local"`. We're only interested in `"air-quality"`, so let's connect to that one:

db = client["air-quality"]

In MongoDB, a **database** is a container for **collections**. Each database gets its own set of files, and a single MongoDB **server** typically has multiple databases.

## Collections

for c in db.list_collections():

system.views
lagos
system.buckets.lagos
nairobi
system.buckets.nairobi
dar-es-salaam
system.buckets.dar-es-salaam

As you can see, there are three actual collections here: `"nairobi"`, `"lagos"`, and `"dar-es-salaam"`. Since we're only interested in the `"lagos"` collection, let's get it on its own like this:

lagos = db["lagos"]

## Documents

A MongoDB **document** is an individual record of data in a **collection**, and is the basic unit of analysis in MongoDB. Documents come with **metadata** that helps us understand what the document is; we'll get back to that in a minute. In the meantime, let's use the [`count_documents`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.count_documents) method to see how many documents the `"lagos"` collection contains:

lagos.count_documents({})

166496

<font size="+1">Practice</font>

client = MongoClient(host = "localhost", port = 27017)

[    {'empty': False, 'name': 'admin', 'sizeOnDisk': 40960},
     {'empty': False, 'name': 'air-quality', 'sizeOnDisk': 7000064},
     {'empty': False, 'name': 'config', 'sizeOnDisk': 98304},
     {'empty': False, 'name': 'local', 'sizeOnDisk': 73728},
     {'empty': False, 'name': 'wqu-abtest', 'sizeOnDisk': 585728}]
system.views
lagos
system.buckets.lagos
nairobi
system.buckets.nairobi
dar-es-salaam
system.buckets.dar-es-salaam

202212

### Retrieving Data

result = lagos.find_one({})

{    '_id': ObjectId('6334b0f18c51459f9b1d955d'),
     'metadata': {    'lat': 6.501,
                      'lon': 3.367,
                      'measurement': 'temperature',
                      'sensor_id': 10,
                      'sensor_type': 'DHT11',
                      'site': 4},
     'temperature': nan,
     'timestamp': datetime.datetime(2018, 1, 7, 7, 7, 3, 88000)}

### Key-Value Pairs

### Metadata

### Identifiers

<font size="+1">Practice</font>

result = nairobi.find_one({})

{    'P1': 39.67,
     '_id': ObjectId('6334b0e98c51459f9b198d27'),
     'metadata': {    'lat': -1.3,
                      'lon': 36.785,
                      'measurement': 'P1',
                      'sensor_id': 57,
                      'sensor_type': 'SDS011',
                      'site': 29},
     'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}

## Analyzing Data

lagos.distinct("metadata.site")

[3, 4]

Notice that in order to grab the `"site"` number, we needed to include the `"metadata"` tag.

print("Documents from site 3:", lagos.count_documents({"metadata.site": 3}))

Documents from site 3: 140586
Documents from site 4: 25910

<font size="+1">Practice</font>

print("Documents from site 29:", nairobi.count_documents({"metadata.site":29}))

Documents from site 29: 131852
Documents from site 6: 70360

print("Documents from site 29:", nairobi.count_documents({"metadata.site": 29}))

Documents from site 29: 131852
Documents from site 6: 70360

Now that we know how many *documents* are associated with each site, let's keep drilling down and find the number of *readings* for each site. We'll do this with the [`aggregate`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.aggregate) method.

    [{"$group": {"_id": "$metadata.site", "count": {"$count": {}}}}]

[{'_id': 4, 'count': 25910}, {'_id': 3, 'count': 140586}]

With that information in mind, we might want to know what those readings actually are. Since we're really interested in measures of air quality, let's take a look at the `P2` values in the `"lagos"` collection. `P2` measures the amount of particulate matter in the air, which in this case is something called PM 2.5. If we wanted to get all the documents in a collection, we could, but that would result in an unmanageably large number of records clogging up the memory on our machines. Instead, let's use the [`find`](https://pymongo.readthedocs.io/en/stable/api/pymongo/collection.html#pymongo.collection.Collection.find) method and use `limit` to make sure we only get back the first 3.

result = lagos.find({"metadata.measurement": "P2"}).limit(3)

[    {    'P2': 14.42,
          '_id': ObjectId('6334b0f28c51459f9b1de145'),
          'metadata': {    'lat': 6.501,
                           'lon': 3.367,
                           'measurement': 'P2',
                           'sensor_id': 6,
                           'sensor_type': 'PPD42NS',
                           'site': 4},
          'timestamp': datetime.datetime(2018, 1, 7, 7, 7, 3, 39000)},
     {    'P2': 19.66,
          '_id': ObjectId('6334b0f28c51459f9b1de146'),
          'metadata': {    'lat': 6.501,
                           'lon': 3.367,
                           'measurement': 'P2',
                           'sensor_id': 6,
                           'sensor_type': 'PPD42NS',
                           'site': 4},
          'timestamp': datetime.datetime(2018, 1, 7, 7, 11, 23, 870000)},
     {    'P2': 24.79,
          '_id': ObjectId('6334b0f28c51459f9b1de147'),
          'metadata': {    'lat': 6.501,
                           'lon': 3.367,
                           'measurement': 'P2',
                           'sensor_id': 6,
                           'sensor_type': 'PPD42NS',
                           'site': 4},
          'timestamp': datetime.datetime(2018, 1, 7, 7, 21, 53, 981000)}]

<font size="+1">Practice</font>

    [{"$group": {"_id": "$metadata.site", "count": {"$count": {}}}}]

[{'_id': 29, 'count': 131852}, {'_id': 6, 'count': 70360}]
[    {    'P2': 14.42,
          '_id': ObjectId('6334b0f28c51459f9b1de145'),
          'metadata': {    'lat': 6.501,
                           'lon': 3.367,
                           'measurement': 'P2',
                           'sensor_id': 6,
                           'sensor_type': 'PPD42NS',
                           'site': 4},
          'timestamp': datetime.datetime(2018, 1, 7, 7, 7, 3, 39000)},
     {    'P2': 19.66,
          '_id': ObjectId('6334b0f28c51459f9b1de146'),
          'metadata': {    'lat': 6.501,
                           'lon': 3.367,
                           'measurement': 'P2',
                           'sensor_id': 6,
                           'sensor_type': 'PPD42NS',
                           'site': 4},
          'timestamp': datetime.datetime(2018, 1, 7, 7, 11, 23, 870000)},
     {    'P2': 24.79,
          '_id': ObjectId('6334b0f28c51459f9b1de147'),
          'metadata': {    'lat': 6.501,
                           'lon': 3.367,
                           'measurement': 'P2',
                           'sensor_id': 6,
                           'sensor_type': 'PPD42NS',
                           'site': 4},
          'timestamp': datetime.datetime(2018, 1, 7, 7, 21, 53, 981000)}]

So far, we've been dealing with relatively small subsets of the data in our collections, but what if we need to work with something bigger? Let's start by using `distinct` to remind ourselves of the kinds of data we have at our disposal.

lagos.distinct("metadata.measurement")

['temperature', 'P1', 'humidity', 'P2']

There are also comparison query operators that can be helpful for filtering the data. In total, we have

result = nairobi.find({"timestamp": {"$gt": datetime.datetime(2018, 9, 1)}}).limit(3)

[    {    'P1': 39.67,
          '_id': ObjectId('6334b0e98c51459f9b198d27'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)},
     {    'P1': 39.13,
          '_id': ObjectId('6334b0e98c51459f9b198d28'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 5, 3, 941000)},
     {    'P1': 30.07,
          '_id': ObjectId('6334b0e98c51459f9b198d29'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 10, 4, 374000)}]

result = nairobi.find({"timestamp": {"$lt": datetime.datetime(2018, 12, 1)}}).limit(3)

[    {    'P1': 39.67,
          '_id': ObjectId('6334b0e98c51459f9b198d27'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)},
     {    'P1': 39.13,
          '_id': ObjectId('6334b0e98c51459f9b198d28'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 5, 3, 941000)},
     {    'P1': 30.07,
          '_id': ObjectId('6334b0e98c51459f9b198d29'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 10, 4, 374000)}]

    {"timestamp": {"$eq": datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}}

[    {    'P1': 39.67,
          '_id': ObjectId('6334b0e98c51459f9b198d27'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)},
     {    'P2': 34.43,
          '_id': ObjectId('6334b0ea8c51459f9b1a0db2'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P2',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}]

<font size="+1">Practice</font>

result = nairobi.find({"timestamp":{"$gte":datetime.datetime(2018,12,1,0,0,6,767000)}}).limit(3)

[    {    'P1': 17.08,
          '_id': ObjectId('6334b0e98c51459f9b19eba8'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 12, 1, 0, 0, 6, 767000)},
     {    'P1': 17.62,
          '_id': ObjectId('6334b0e98c51459f9b19eba9'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 12, 1, 0, 5, 6, 327000)},
     {    'P1': 11.05,
          '_id': ObjectId('6334b0e98c51459f9b19ebaa'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 12, 1, 0, 10, 5, 579000)}]

# Less than or equal to

## Updating Documents

We can also update documents by passing some filter and new values using `update_one` to update one record or `update_many` to update many records. Let's look at an example. Before updating, we have this record showing like this:

    {"timestamp": {"$eq": datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}}

[    {    'P1': 39.67,
          '_id': ObjectId('6334b0e98c51459f9b198d27'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}]

Now we are updating the sensor type from `"SDS011"` to `"SDS"`, we first select all records with sensor type equal to `"SDS011"`, then set the new value to `"SDS"`:

    {"metadata.sensor_type": {"$eq": "SDS101"}},

Now we can see all records have changed:

    {"timestamp": {"$eq": datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}}

[    {    'P1': 39.67,
          '_id': ObjectId('6334b0e98c51459f9b198d27'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P1',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)},
     {    'P2': 34.43,
          '_id': ObjectId('6334b0ea8c51459f9b1a0db2'),
          'metadata': {    'lat': -1.3,
                           'lon': 36.785,
                           'measurement': 'P2',
                           'sensor_id': 57,
                           'sensor_type': 'SDS011',
                           'site': 29},
          'timestamp': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}]

We can change it back:

    {"$set": {"metadata.sensor_type": "SDS101"}},

result.raw_result

{'n': 0, 'nModified': 0, 'ok': 1.0, 'updatedExisting': False}

## Aggregation

Since we're looking for *big* numbers, we need to figure out which one of those dimensions has the largest number of measurements by **aggregating** the data in each document. Since we already know that `site 3` has significantly more documents than `site 2`, let's start looking at `site 3`. We can use the `$match` syntax to only select `site 3` data. The code to do that looks like this:

        {"$group": {"_id": "$metadata.measurement", "count": {"$count": {}}}},

[    {'_id': 'temperature', 'count': 35147},
     {'_id': 'P1', 'count': 35146},
     {'_id': 'P2', 'count': 35146},
     {'_id': 'humidity', 'count': 35147}]

<font size="+1">Practice</font>

pp.pprint(list(result))

After aggregation, there is another useful operator called `$project`, which allows you to specify which fields to display by adding new fields or deleting fields. Using the Nairobi data from site 29, we can first count each sensor type:<span style='color: transparent; font-size:1%'>WQU WorldQuant University Applied Data Science Lab QQQQ</span>

        {"$group": {"_id": "$metadata.sensor_type", "count": {"$count": {}}}},

[{'_id': 'DHT22', 'count': 66038}, {'_id': 'SDS011', 'count': 65814}]

We can see there are two sensor types and the corresponding counts. If we only want to display what are the types but do not care about the counts, we can suppress the `count` filed by setting it at 0 in `$project`:

        {"$group": {"_id": "$metadata.sensor_type", "count": {"$count": {}}}},

[{'_id': 'DHT22'}, {'_id': 'SDS011'}]

The `$project` syntax is also useful for deleting the intermediate fields that we used to generate our final fields but no longer need. In the following example, let's calculate the date difference for each sensor type. We'll first use the aggregation method to get the start date and last date.

                "date_min": {"$min": "$timestamp"},

[    {    '_id': 'DHT22',
          'date_max': datetime.datetime(2018, 12, 31, 23, 57, 37, 257000),
          'date_min': datetime.datetime(2018, 9, 1, 0, 0, 4, 301000)},
     {    '_id': 'SDS011',
          'date_max': datetime.datetime(2018, 12, 31, 23, 55, 4, 811000),
          'date_min': datetime.datetime(2018, 9, 1, 0, 0, 2, 472000)}]

Then we can calculate the date difference using `$dateDiff`, which gets the date difference through specifying the start date, end date and unit for timestamp data. We can see from the results above that the dates, are very close to each other. The only differences are in the minutes, so we can specify the unit as minute to show the difference. Since we don't need the start date and end dates, we can define a `"dateDiff"` field inside `$project`, so that it will be shown in the final display:

                "date_min": {"$min": "$timestamp"},

[{'_id': 'DHT22', 'dateDiff': 175677}, {'_id': 'SDS011', 'dateDiff': 175675}]

If we specify unit as `day`, it will show the difference between the dates:

                "date_min": {"$min": "$timestamp"},

<font size="+1">Practice</font>

pp.pprint(list(result))

We can do more with the date data using `$dateTrunc`, which truncates datetime data. We need to specify the datetime data, which can be a `Date`, a `Timestamp`, or an `ObjectID`. Then we need to specify the `unit` (year, month, day, hour, minute, second) and `binSize` (numerical variable defining the size of the truncation). Let's check the example below, where we group data by the month using `$dateTrunc` and then count how many observations there are for each month.

                            "date": "$timestamp",

<font size="+1">Practice</font>

pp.pprint(list(result))

## Finishing Up

    # `projection` limits the kinds of data we'll get back.

{'P2': 0.01, 'timestamp': datetime.datetime(2018, 10, 1, 0, 0, 52, 906000)}

Finally, we'll use pandas to read the extracted data into a DataFrame, making sure to set `timestamp` as the index:

df = pd.DataFrame(result).set_index("timestamp")

<font size="+1">Practice</font>

result = ...

# References & Further Reading

---

​x

<font size="+3"><strong>3.4. ARMA Models</strong></font>

ximport inspectimport timeimport warnings​import matplotlib.pyplot as pltimport pandas as pdimport plotly.express as pximport seaborn as snsfrom IPython.display import VimeoVideofrom pymongo import MongoClientfrom sklearn.metrics import mean_absolute_errorfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacffrom statsmodels.tsa.arima.model import ARIMA​warnings.filterwarnings("ignore")

/opt/conda/lib/python3.9/site-packages/statsmodels/tsa/base/tsa_model.py:7: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  from pandas import (to_datetime, Int64Index, DatetimeIndex, Period,
/opt/conda/lib/python3.9/site-packages/statsmodels/tsa/base/tsa_model.py:7: FutureWarning: pandas.Float64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  from pandas import (to_datetime, Int64Index, DatetimeIndex, Period,

 VimeoVideo("665851728", h="95c59d2805", width=600)

# Prepare Data

## Import

**Task 3.4.1:** Create a client to connect to the MongoDB server, then assign the `"air-quality"` database to `db`, and the `"nairobi"` collection to `nairobi`.

 client = MongoClient(host="localhost",port=27017)db = client["air-quality"]nairobi = db["nairobi"]

 def wrangle(collection,resample_rule="1H"):​    results = collection.find(        {"metadata.site": 29, "metadata.measurement": "P2"},        projection={"P2": 1, "timestamp": 1, "_id": 0},    )​    # Read results into DataFrame    df = pd.DataFrame(list(results)).set_index("timestamp")​    # Localize timezone    df.index = df.index.tz_localize("UTC").tz_convert("Africa/Nairobi")​    # Remove outliers    df = df[df["P2"] < 500]​    # Resample and forward-fill    y = df["P2"].resample(resample_rule).mean().fillna(method="ffill")​    return y

 VimeoVideo("665851670", h="3efc0c20d4", width=600)

**Task 3.4.2:** Change your `wrangle` function so that it has a `resample_rule` argument that allows the user to change the resampling interval. The argument default should be `"1H"`.

 # Check your workfunc_params = set(inspect.signature(wrangle).parameters.keys())assert func_params == set(    ["collection", "resample_rule"]), f"Your function should take two arguments: `'collection'`, `'resample_rule'`. Your function takes the following arguments: {func_params}"

**Task 3.4.3:** Use your wrangle function to read the data from the `nairobi` collection into the Series `y`.

 y = wrangle(nairobi,resample_rule="1H")y.head()

timestamp
2018-09-01 03:00:00+03:00    17.541667
2018-09-01 04:00:00+03:00    15.800000
2018-09-01 05:00:00+03:00    11.420000
2018-09-01 06:00:00+03:00    11.614167
2018-09-01 07:00:00+03:00    17.665000
Freq: H, Name: P2, dtype: float64

 # Check your workassert isinstance(y, pd.Series), f"`y` should be a Series, not a {type(y)}."assert len(y) == 2928, f"`y` should have 2,928 observations, not {len(y)}."assert (    y.isnull().sum() == 0), f"There should be no null values in `y`. Your `y` has {y.isnull().sum()} null values."

## Explore

 VimeoVideo("665851654", h="687ff8d5ee", width=600)

**Task 3.4.4:** Create an ACF plot for the data in `y`. Be sure to label the x-axis as `"Lag [hours]"` and the y-axis as `"Correlation Coefficient"`.

 fig,ax=plt.subplots(figsize=(15,6))plot_acf(y,ax=ax)plt.xlabel("lag in hours")plt.ylabel("Correlation coefficient")

Text(0, 0.5, 'Correlation coefficient')

 VimeoVideo("665851644", h="e857f05bfb", width=600)

**Task 3.4.5:** Create an PACF plot for the data in `y`. Be sure to label the x-axis as `"Lag [hours]"` and the y-axis as `"Correlation Coefficient"`.

 fig,ax = plt.subplots(figsize=(15,6))plot_pacf(y,ax=ax)plt.xlabel("Lag [hours]")plt.ylabel("Correlation coefficient")

Text(0, 0.5, 'Correlation coefficient')

## Split

**Task 3.4.6:** Create a training set `y_train` that contains only readings from October 2018, and a test set `y_test` that contains readings from November 1, 2018.

 # print(len(y))y_train = y.get("October 2018")y_test = y.get("November 1, 2018")

 print(len(y_test))

24

 # Check your workassert (    len(y_train) == 744), f"`y_train` should have 744 observations, not {len(y_train)}."assert len(y_test) == 24, f"`y_test` should have 24 observations, not {len(y_test)}."

# Build Model

## Baseline

**Task 3.4.7:** Calculate the baseline mean absolute error for your model.

 y_train_mean = y_train.mean()y_pred_baseline = [y_train_mean] * len(y_train)mae_baseline = mean_absolute_error(y_train, y_pred_baseline)print("Mean P2 Reading:", round(y_train_mean, 2))print("Baseline MAE:", round(mae_baseline, 2))

Mean P2 Reading: 10.12
Baseline MAE: 4.17

## Iterate

 VimeoVideo("665851576", h="36e2dc6269", width=600)

**Task 3.4.8:** Create ranges for possible $p$ and $q$ values. `p_params` should range between 0 and 25, by steps of 8. `q_params` should range between 0 and 3 by steps of 1.

 p_params = range(0,25,8)q_params = range(0,3,1)list(p_params)

[0, 8, 16, 24]

 VimeoVideo("665851476", h="d60346ed30", width=600)

**Task 3.4.9:** Complete the code below to train a model with every combination of hyperparameters in `p_params` and `q_params`. Every time the model is trained, the mean absolute error is calculated and then saved to a dictionary. If you're not sure where to start, do the code-along with Nicholas!

 # Create dictionary to store MAEsmae_grid = dict()# Outer loop: Iterate through possible values for `p`for p in p_params:    # Create key-value pair in dict. Key is `p`, value is empty list.    mae_grid[p] = list()    # Inner loop: Iterate through possible values for `q`    for q in q_params:        # Combination of hyperparameters for model        order = (p, 0, q)        # Note start time        start_time = time.time()        # Train model        model = ARIMA(y_train, order=order).fit()        # Calculate model training time        elapsed_time = round(time.time() - start_time, 2)        print(f"Trained ARIMA {order} in {elapsed_time} seconds.")        # Generate in-sample (training) predictions        y_pred = model.predict()        # Calculate training MAE        mae = mean_absolute_error(y_train,y_pred)        # Append MAE to list in dictionary        mae_grid[p].append(mae)​print()print(mae_grid)

Trained ARIMA (0, 0, 0) in 0.55 seconds.
Trained ARIMA (0, 0, 1) in 0.35 seconds.
Trained ARIMA (0, 0, 2) in 1.1 seconds.
Trained ARIMA (8, 0, 0) in 10.3 seconds.
Trained ARIMA (8, 0, 1) in 42.8 seconds.
Trained ARIMA (8, 0, 2) in 80.41 seconds.
Trained ARIMA (16, 0, 0) in 42.9 seconds.
Trained ARIMA (16, 0, 1) in 159.28 seconds.
Trained ARIMA (16, 0, 2) in 220.81 seconds.
Trained ARIMA (24, 0, 0) in 122.58 seconds.
Trained ARIMA (24, 0, 1) in 158.0 seconds.
Trained ARIMA (24, 0, 2) in 288.89 seconds.

{0: [4.171460443827197, 3.3506427433555537, 3.1057222587888647], 8: [2.9384480570404223, 2.9149008203151663, 2.908046094677133], 16: [2.9201084726122, 2.929436207635203, 2.913638894139686], 24: [2.914390327277196, 2.9136013252035013, 2.89792511429827]}

 VimeoVideo("665851464", h="12f4080d0b", width=600)

**Task 3.4.10:** Organize all the MAE's from above in a DataFrame names `mae_df`. Each row represents a possible value for $q$ and each column represents a possible value for $p$.

 mae_gridmae_df = pd.DataFrame(mae_grid)mae_df.round(4)

 VimeoVideo("665851453", h="dfd415bc08", width=600)

**Task 3.4.11:** Create heatmap of the values in `mae_grid`. Be sure to label your x-axis `"p values"` and your y-axis `"q values"`.

 sns.heatmap(mae_df,cmap="Blues")plt.xlabel("p values")plt.ylabel("q values")​

Text(33.0, 0.5, 'q values')

 VimeoVideo("665851444", h="8b58161f26", width=600)

**Task 3.4.12:** Use the [`plot_diagnostics`](https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima.model.ARIMAResults.plot_diagnostics.html) method to check the residuals for your model. Keep in mind that the plot will represent the residuals from the last model you trained, so make sure it was your best model, too!

 fig, ax = plt.subplots(figsize=(15, 12))model.plot_diagnostics(fig=fig)​

## Evaluate

 VimeoVideo("665851439", h="c48d80cdf4", width=600)

**Task 3.4.13:** Complete the code below to perform walk-forward validation for your model for the entire test set `y_test`. Store your model's predictions in the Series `y_pred_wfv`. Choose the values for $p$ and $q$ that best balance model performance and computation time. Remember: This model is going to have to train 24 times before you can see your test MAE!<span style='color: transparent; font-size:1%'>WQU WorldQuant University Applied Data Science Lab QQQQ</span>

 y_pred_wfv = pd.Series()history = y_train.copy()for i in range(len(y_test)):    model = ARIMA(history,order=(8,0,1)).fit()    next_pred = model.forecast()    y_pred_wfv = y_pred_wfv.append(next_pred)    history = history.append(y_test[next_pred.index])y_pred_wfv

2018-11-01 00:00:00+03:00     8.725735
2018-11-01 01:00:00+03:00     7.988511
2018-11-01 02:00:00+03:00     7.133804
2018-11-01 03:00:00+03:00     6.372284
2018-11-01 04:00:00+03:00     7.539867
2018-11-01 05:00:00+03:00     8.538572
2018-11-01 06:00:00+03:00    10.757824
2018-11-01 07:00:00+03:00    10.611658
2018-11-01 08:00:00+03:00     9.914120
2018-11-01 09:00:00+03:00    10.911411
2018-11-01 10:00:00+03:00     9.094045
2018-11-01 11:00:00+03:00     7.250349
2018-11-01 12:00:00+03:00     6.515896
2018-11-01 13:00:00+03:00     7.184177
2018-11-01 14:00:00+03:00     6.864486
2018-11-01 15:00:00+03:00     7.564168
2018-11-01 16:00:00+03:00     8.732164
2018-11-01 17:00:00+03:00     8.295855
2018-11-01 18:00:00+03:00     8.199158
2018-11-01 19:00:00+03:00     8.293687
2018-11-01 20:00:00+03:00     9.204907
2018-11-01 21:00:00+03:00     9.629578
2018-11-01 22:00:00+03:00     9.658772
2018-11-01 23:00:00+03:00    10.351504
Freq: H, dtype: float64

  test_mae = mean_absolute_error(y_test,y_pred_wfv)print("Test MAE (walk forward validation):", round(test_mae, 2))

Test MAE (walk forward validation): 1.67

# Communicate Results

 VimeoVideo("665851423", h="8236ff348f", width=600)

**Task 3.4.14:** First, generate the list of training predictions for your model. Next, create a DataFrame `df_predictions` with the true values `y_test` and your predictions `y_pred_wfv` (don't forget the index). Finally, plot `df_predictions` using plotly express. Make sure that the y-axis is labeled `"P2"`.

 df_predictions = pd.DataFrame(    {"y_test":y_test,"y_pred_wfv":y_pred_wfv}​)fig = px.line(df_predictions,labels={"Value":"PM2.5"})fig.show()

---

xxxxxxxxxx​

xxxxxxxxxx<font size="+3"><strong>Pandas: Advanced</strong></font>

xxxxxxxxxx# Calculate Summary Statistics for a DataFrame or Series

xxxxxxxxxxMany datasets are large, and it can be helpful to get a summary of information in them. Let's load a dataset and examine the first few rows:

mexico_city1 = pd.read_csv("./data/mexico-city-real-estate-1.csv")

xxxxxxxxxxLet's get a summary description of this dataset.

mexico_city1.describe()

xxxxxxxxxxLike most large datasets, this one has many values which are missing. The describe function will ignore missing values in each column. You can also remove rows and columns with missing values, and then get a summary of the data that's still there. We need to remove columns first, before removing the rows; the sequence of operations here is important. The code looks like this:

    ["floor", "price_usd_per_m2", "expenses", "rooms"], axis=1

xxxxxxxxxxLet's take a look at our new, cleaner dataset.

mexico_city1.describe()

xxxxxxxxxx<font size="+1">Practice</font>

mexico_city1 = pd.read_csv("./data/mexico-city-real-estate-1.csv")

xxxxxxxxxx# Select a Series from a DataFrame

xxxxxxxxxxSince the datasets we work with are so large, you might want to focus on a single column of a DataFrame. Let's load up the `mexico-city-real-estate-2` dataset, and examine the first few rows to find the column names.

mexico_city2 = pd.read_csv("./data/mexico-city-real-estate-2.csv")

xxxxxxxxxxMaybe we're interested in the `surface_covered_in_m2` column. The code to extract just that one column looks like this:

surface_covered_in_m2 = mexico_city2["surface_covered_in_m2"]

xxxxxxxxxx<font size="+1">Practice</font>

print(price)

xxxxxxxxxx# Subset a DataFrame by Selecting One or More Columns

xxxxxxxxxxYou may find it more efficient to work with a smaller portion of a dataset that's relevant to you. One way to do this is to select some columns from a DataFrame and make a new DataFrame with them. Let's load a dataset to do this and examine the first few rows to find the column headings:

mexico_city4 = pd.read_csv("./data/mexico-city-real-estate-4.csv")

xxxxxxxxxxLet's choose `"operation"`, `"property_type"`, `"place_with_parent_names"`, and `"price"`:

    ["operation", "property_type", "place_with_parent_names", "price"]

xxxxxxxxxxOnce we've done that, we can find the resulting number of entries in the DataFrame:

mexico_city4_subset.shape

xxxxxxxxxx<font size="+1">Practice</font>

    ["operation", "lat-lon", "place_with_parent_names"]

xxxxxxxxxx# Subset the Columns of a DataFrame Based on Data Types

xxxxxxxxxxIt's helpful to be able to find specific types of entries — typically numeric ones — and put all of these in a separate DataFrame. First, let's take a look at the `mexico-city-real-estate-5` dataset.

mexico_city5 = pd.read_csv("./data/mexico-city-real-estate-5.csv")

xxxxxxxxxxThe code to subset just the numerical entries looks like this:

mexico_city5_numbers = mexico_city5.select_dtypes(include="number")

xxxxxxxxxx<font size="+1">Practice</font>

print(mexico_city3_no_numbers.shape)

xxxxxxxxxx# Working with `value_counts` in a Series

df1 = pd.read_csv("data/colombia-real-estate-1.csv", usecols=["department"])

xxxxxxxxxx<font size="+1">Practice</font>

df2 = ...

xxxxxxxxxx# Series and `Groupby`

dept_group = df1.groupby("department")

xxxxxxxxxxTo make sure we got all the departments in the dataset, let's print the first occurrence of each department.

dept_group.first()

xxxxxxxxxx<font size="+1">Practice</font>

dept_group.first()

xxxxxxxxxxNow that we have all the properties grouped by department, we might want to see the properties in just one of the departments. We can use the `get_group` method to do that. If we just wanted to see the properties in `"Santander"`, for example, the code would look like this:

dept_group1 = df1.groupby("department")

xxxxxxxxxxWe can also make groups based on more than one category by adding them to the `groupby` method. After resetting the `df1` DataFrame, here's what the code looks like if we want to group properties both by `department` and by `property_type`.

dept_group2 = df1.groupby(["department", "property_type"])

xxxxxxxxxx<font size="+1">Practice</font>

dept_group.first()

xxxxxxxxxxFinally, it's possible to use `groupby` to calculate aggregations. For example, if we wanted to find the average property area in each department, we would use the `.mean()` method. This is what the code for that looks like:

dept_group = df1.groupby("department")["area_m2"].mean()

xxxxxxxxxx*Practice*

dept_group = ...

xxxxxxxxxx# Subset a DataFrame's Columns Based on the Column Data Types

xxxxxxxxxxIt's helpful to be able to find entries of a certain type, typically numerical entries, and put all of these in a separate DataFrame. Let's load a dataset to see how this works:

mexico_city5 = pd.read_csv("./data/mexico-city-real-estate-5.csv")

xxxxxxxxxxNow let's get only numerical entries:

mexico_city5_numbers = mexico_city5.select_dtypes(include="number")

xxxxxxxxxxLet's now find all entries which are not numerical entries:

mexico_city5_no_numbers = mexico_city5.select_dtypes(exclude="number")

xxxxxxxxxx<font size="+1">Practice</font>

print(mexico_city3_no_numbers.shape)

xxxxxxxxxx# Subset a DataFrame's columns based on column names

xxxxxxxxxxTo subset a DataFrame by column names, either define a list of columns to include or define a list of columns to exclude. Once you've done that, you can retain or drop the columns accordingly. For example, let's suppose we want to modify the `mexico_city3` dataset and only retain the first three columns. Let's define two lists, one with the columns to retain and one with the columns to drop:

keep_cols = ["operation", "property_type", "place_with_parent_names"]

xxxxxxxxxxNext, we'll explore both approaches to subset `mexico_city3`. To retain columns based on `keep_cols`:

mexico_city3_subsetted = mexico_city3[keep_cols]

xxxxxxxxxxTo drop columns in `drop_cols`:

mexico_city3_subsetted = mexico_city3.drop(columns=drop_cols)

xxxxxxxxxx<font size="+1">Practice</font>

​

xxxxxxxxxx## Pivot Tables

mexico_city3_pivot = ...

xxxxxxxxxx# Subsetting with Masks

df1 = pd.read_csv("data/colombia-real-estate-1.csv")

xxxxxxxxxxNotice that `mask` is a Series of Boolean values. Where properties are smaller than 200 square meters, our statement evaluates as `False`; where they're bigger than 200, it evaluates to `True`.

df1[mask].head()

xxxxxxxxxx<font size="+1">Practice</font>

df2[mask].head()

xxxxxxxxxxWe can also create masks with multiple criteria using `&` for "and" and `|` for "or." For example, here's how we would find all properties in Atlántico that are bigger than 400 square meters.

mask = (df1["area_m2"] > 400) & (df1["department"] == "Atlántico")

xxxxxxxxxx<font size="+1">Practice</font>

df2[mask].head()

xxxxxxxxxx## Reshape a DataFrame based on column values

xxxxxxxxxx## What's a pivot table?

pivot1 = pd.pivot_table(df2, values="price_cop", index="property_type", aggfunc=np.mean)

xxxxxxxxxx<font size="+1">Practice</font>

pivot2 = pd.pivot_table(df2, values="price_cop", index="department", aggfunc=np.mean)

xxxxxxxxxx## Combine multiple categories in a Series

xxxxxxxxxxCategorical variables can be collapsed into a fewer number of categories. One approach is to retain the values of the most frequently observed values and collapse all remaining values into a single category. For example, to retain only the values of the top 10 most frequent categories in the `department` column and then collapse the other categories together, use `value_counts` to generate the count of the values:

df2["department"].value_counts()

xxxxxxxxxxNext, select just the top 10 using the `head()` method, and select the column names by using the `index` attribute of the series:

top_10 = df2["department"].value_counts().head(10).index

xxxxxxxxxxFinally, use the apply method and a lambda function to select only the values from the `department` column and collapse the remaining values into the value `Other`:

df2["department"] = df2["department"].apply(lambda x: x if x in top_10 else "Other")

xxxxxxxxxx<font size="+1">Practice</font>

​

xxxxxxxxxx# Cross Tabulation

xxxxxxxxxxThe pandas `crosstab` function is a useful for working with grouped summary statistics for **categorical** data. It starts by picking two categorical columns, then defines one as the index and the other as the column. If the aggregate function and value column is not defined, `crosstab` will simply calculate the frequency of each combination by default. Let's see the example below from the Colombia real estate dataset.

df = pd.read_csv("data/colombia-real-estate-1.csv")

xxxxxxxxxxThe following function calculates the frequency of each combination for two variables, `department` and `property_type`, in the data set.

pd.crosstab(index=df["department"], columns=df["property_type"])

xxxxxxxxxxFrom the previous example, you can see we have created a DataFrame with the index showing unique observation in variable `department`, while the columns are the unique observation for the variable `property_type`. Each cell shows how many data points there are for each combination of `department` type and `property_type`. For example, there are 8 `apartments` in department `Antioquia`.

xxxxxxxxxxWe can further specify a value column and an aggregate function, like in `pivot_table`, to conduct more complicated calculations for the two categorical variables. In the following example, we're looking at the average area size for different property types in the departments in Colombia.

    columns=df["property_type"],

xxxxxxxxxx<font size="+1">Practice</font>

# Create `crosstab`

xxxxxxxxxx# Applying Functions to DataFrames and Series

xxxxxxxxxx`apply` is a useful method for to using one function on all the rows or columns of a DataFrame efficiently. Let's take the following real estate dataset as an example:

df = pd.read_csv("data/colombia-real-estate-2.csv", usecols=["area_m2", "price_cop"])

xxxxxxxxxxBy specifying the function inside `apply()`, we can transform the whole DataFrame. For example, I am calculating the square root of each row at each column:

import numpy as np

xxxxxxxxxxNote you can also specify the `"axis"` inside `apply`. By default, we have `axis=0`, which means we are applying the function to each column. We can also switch to `axis=1` if we want to apply the function to each row. See the following example showing the sum of all rows for each column:

df.apply(np.sum, axis=0)

xxxxxxxxxxThe following code will get the sum of all columns for each row:

df.apply(np.sum, axis=1)

xxxxxxxxxxBy specifying the column name, we can also apply the function to a specific column or columns. Note that we can also specify index (row names) to only apply functions to specific rows, however, it is not common in practice.

df["area_m2"].apply(np.sqrt)

xxxxxxxxxxWe can assign the results to a new column:

df["area_sqrt"] = df["area_m2"].apply(np.sqrt)

xxxxxxxxxx<font size="+1">Practice</font>

df["sum_columns"] = ...

xxxxxxxxxx`df.aggregate()`, or `df.agg()`, shares the same concept as `df.apply()` in terms of applying functions to a DataFrame, but `df.aggregate()` can only apply aggregate functions like `sum`, `mean`, `min`, `max`, etc. See the following example for more details:

df = pd.read_csv("data/colombia-real-estate-2.csv", usecols=["area_m2", "price_cop"])

xxxxxxxxxxWe can check what's the minimum number for each column calling the `min` aggregate function:

df.agg("min")

xxxxxxxxxxLike `apply()`, we can also specify the `axis` argument to switch axis:

df.agg("min", axis=1)

xxxxxxxxxxWe can apply aggregate function to a whole DataFrame using `df.agg()`, or specify the column name for a subset of DataFrame:

df["area_m2"].agg("min")

xxxxxxxxxxWe can also apply a list of aggregate functions to a DataFrame:

df.agg(["sum", "max"])

xxxxxxxxxxThe syntax above will calculate both `sum` and `max` for each column, and store the result as index. Besides, we can also apply different aggregate functions to different columns. In this case, we need to pass a dictionary specifying key as column names, and value as corresponding aggregate function names:

df.agg({"area_m2": "sum", "price_cop": "min"})

xxxxxxxxxx<font size="+1">Practice</font>

​

xxxxxxxxxx# Working with Dates

xxxxxxxxxx## Time Stamps

xxxxxxxxxxPandas' time series capabilities are built on the `Timestamp` class. For instance, we can transform date strings of different formats into `Timestamp`:

print(pd.Timestamp("January 8, 2022"))

xxxxxxxxxxWe can also do date math with `Timestamp`s. Note that when we subtract two dates, we get a special object called a `Timedelta`.

time_delta = pd.Timestamp("Feb. 11 2016 2:30 am") - pd.Timestamp("2015-08-03 5:14 pm")

xxxxxxxxxxThe pandas `Timestamp` class is also compatible with Python's `datetime` module.<span style='color: transparent; font-size:1%'>WQU WorldQuant University Applied Data Science Lab QQQQ</span>

pd.Timestamp("Feb 11, 2016") - datetime.datetime(2015, 8, 3)

xxxxxxxxxxWe can manipulate dates using some function from `offset()`. We can use the `Day()` function to add or decrease days from a `timestamp`. The following function calculates the date that's four days prior to 9 January 2017:

from pandas.tseries.offsets import BDay, BMonthEnd, Day

xxxxxxxxxxIn some case, you might need to add or subtract only business days. That's when you use `BDay()`:

print(pd.Timestamp("January 9, 2017") - BDay(4))

xxxxxxxxxxWe can also do an offset at the monthly level. For example, you can use `BMonthEnd(n)` to find the last business day $n$ months later. The following function shows the last business day of January 2017:

print(pd.Timestamp("January 9, 2017") + BMonthEnd(0))

xxxxxxxxxxThe following function shows the last business day for May 2017, which is four months after February:

print(pd.Timestamp("February 9, 2017") + BMonthEnd(4))

xxxxxxxxxxWe can also use the `strptime` function inside `datetime` to transform string to date format:

date = datetime.strptime("16 Oct 2022 12:14:05", "%d %b %Y %H:%M:%S")

xxxxxxxxxxWe can further transform it into the [ISO 8601 format](https://en.wikipedia.org/wiki/ISO_8601):

date.isoformat()

xxxxxxxxxx## Date Time Indices

xxxxxxxxxx`DatetimeIndex` is a convenient function that transforms date-like data into readable `Datetime` format for a DataFrame index. That way, we can better plot and model time series data. Let's check an example. Here we have Apple's stock prices from 1999 10 2022.

data.set_index("date", inplace=True)

xxxxxxxxxxEven though the index looks like *dates*, it is not in *date format*. So when we plot the data, the index doesn't follow the sequence of years, and the x-axis ticks are hard to read.

data["open"].plot()

xxxxxxxxxxYou can see the index is not in date format.

data.index[:5]

xxxxxxxxxxWe can use the `DatetimeIndex` function to transform it into date:

pd.DatetimeIndex(data.index)

xxxxxxxxxxAnd we can set it as the index:

data.index = pd.DatetimeIndex(data.index)

xxxxxxxxxxNow we can see the benefit from plotting:

data["open"].plot()

xxxxxxxxxx## Date Ranges

xxxxxxxxxxIf we're entering time series data into a DataFrame, it will often be useful to create a range of dates using `date_range`. We can create it with different frequencies by specifying `freq`. Here are the days in a specific range:

pd.date_range(start="1/8/2022", end="3/2/2022", freq="D")

xxxxxxxxxxHere are the business dates for a specific range:

pd.date_range(start="1/8/2022", end="3/2/2022", freq="B")

xxxxxxxxxx# References & Further Reading

xxxxxxxxxx---

xxxxxxxxxx​

xxxxxxxxxx<font size="+3"><strong>3.5. Air Quality in Dar es Salaam 🇹🇿</strong></font>

warnings.simplefilter(action="ignore", category=FutureWarning)

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

pp = PrettyPrinter(indent = 2)

xxxxxxxxxx# Prepare Data

xxxxxxxxxx## Connect

xxxxxxxxxx**Task 3.5.1:** Connect to MongoDB server running at host `"localhost"` on port `27017`. Then connect to the `"air-quality"` database and assign the collection for Dar es Salaam to the variable name `dar`.

client = MongoClient(host="localhost",port=27017)

system.views
lagos
system.buckets.lagos
nairobi
system.buckets.nairobi
dar-es-salaam
system.buckets.dar-es-salaam

wqet_grader.grade("Project 3 Assessment", "Task 3.5.1", [dar.name])

xxxxxxxxxx## Explore

xxxxxxxxxx**Task 3.5.2:** Determine the numbers assigned to all the sensor sites in the Dar es Salaam collection. Your submission should be a list of integers. <span style='color: transparent; font-size:1%'>WQU WorldQuant University Applied Data Science Lab QQQQ</span>

sites = dar.distinct("metadata.site")

[23, 11]

wqet_grader.grade("Project 3 Assessment", "Task 3.5.2", sites)

xxxxxxxxxx**Task 3.5.3:** Determine which site in the Dar es Salaam collection has the most sensor readings (of any type, not just PM2.5 readings). You submission `readings_per_site` should be a list of dictionaries that follows this format:

        {"$group":{"_id":"$metadata.site","count":{"$count":{}}}}

[{'_id': 23, 'count': 60020}, {'_id': 11, 'count': 138412}]

wqet_grader.grade("Project 3 Assessment", "Task 3.5.3", readings_per_site)

xxxxxxxxxx## Import

xxxxxxxxxx**Task 3.5.4:** (5 points) Create a `wrangle` function that will extract the PM2.5 readings from the site that has the most total readings in the Dar es Salaam collection. Your function should do the following steps:

result = dar.find_one({})

{ 'P1': 24.3,
  '_id': ObjectId('6334b0ef8c51459f9b1bffb7'),
  'metadata': { 'lat': -6.818,
                'lon': 39.285,
                'measurement': 'P1',
                'sensor_id': 29,
                'sensor_type': 'SDS011',
                'site': 11},
  'timestamp': datetime.datetime(2018, 1, 1, 0, 0, 4, 53000)}

    df.index = df.index.tz_localize("UTC").tz_convert("Africa/Dar_es_Salaam")

xxxxxxxxxxUse your `wrangle` function to query the `dar` collection and return your cleaned results.

y = wrangle(dar)

timestamp
2018-01-01 03:00:00+03:00    9.456327
2018-01-01 04:00:00+03:00    9.400833
2018-01-01 05:00:00+03:00    9.331458
2018-01-01 06:00:00+03:00    9.528776
2018-01-01 07:00:00+03:00    8.861250
Freq: H, Name: P2, dtype: float64

wqet_grader.grade("Project 3 Assessment", "Task 3.5.4", wrangle(dar))

xxxxxxxxxx## Explore Some More

xxxxxxxxxx**Task 3.5.5:** Create a time series plot of the readings in `y`. Label your x-axis `"Date"` and your y-axis `"PM2.5 Level"`. Use the title `"Dar es Salaam PM2.5 Levels"`.

y.plot(xlabel="Date",ylabel="PM2.5 Level",title="Dar es Salaam PM2.5 Levels",ax=ax)

    wqet_grader.grade("Project 3 Assessment", "Task 3.5.5", file)

xxxxxxxxxx**Task 3.5.6:** Plot the rolling average of the readings in `y`. Use a window size of `168` (the number of hours in a week). Label your x-axis `"Date"` and your y-axis `"PM2.5 Level"`. Use the title `"Dar es Salaam PM2.5 Levels, 7-Day Rolling Average"`.

y.rolling(168).mean().plot(xlabel = "Date",ylabel="PM2.5 Level",ax=ax);

    wqet_grader.grade("Project 3 Assessment", "Task 3.5.6", file)

xxxxxxxxxx**Task 3.5.7:** Create an ACF plot for the data in `y`. Be sure to label the x-axis as `"Lag [hours]"` and the y-axis as `"Correlation Coefficient"`. Use the title `"Dar es Salaam PM2.5 Readings, ACF"`.

plt.title("Dar es Salaam PM2.6 Readings, ACF")

    wqet_grader.grade("Project 3 Assessment", "Task 3.5.7", file)

xxxxxxxxxx**Task 3.5.8:** Create an PACF plot for the data in `y`. Be sure to label the x-axis as `"Lag [hours]"` and the y-axis as `"Correlation Coefficient"`. Use the title `"Dar es Salaam PM2.5 Readings, PACF"`.

plt.title("Dar es Salaam PM2.6 Readings, ACF")

    wqet_grader.grade("Project 3 Assessment", "Task 3.5.8", file)

xxxxxxxxxx## Split

xxxxxxxxxx**Task 3.5.9:** Split `y` into training and test sets. The first 90% of the data should be in your training set. The remaining 10% should be in the test set.

print("y_train shape:", y_train.shape)

y_train shape: (1533,)
y_test shape: (171,)

wqet_grader.grade("Project 3 Assessment", "Task 3.5.9a", y_train)

wqet_grader.grade("Project 3 Assessment", "Task 3.5.9b", y_test)

xxxxxxxxxx# Build Model

xxxxxxxxxx## Baseline

xxxxxxxxxx**Task 3.5.10:** Establish the baseline mean absolute error for your model.

mae_baseline = mean_absolute_error(y_train,y_pred_baseline)

Mean P2 Reading: 8.617582545265428
Baseline MAE: 4.07658759405218

wqet_grader.grade("Project 3 Assessment", "Task 3.5.10", [mae_baseline])

xxxxxxxxxx## Iterate

xxxxxxxxxx**Task 3.5.11:** You're going to use an AR model to predict PM2.5 readings, but which hyperparameter settings will give you the best performance? Use a `for` loop to train your AR model on using settings for `p` from 1 to 30. Each time you train a new model, calculate its mean absolute error and append the result to the list `maes`. Then store your results in the Series `mae_series`.

        mae = mean_absolute_error(y_train.iloc[p:],y_pred)

1    0.947888
2    0.933894
3    0.920850
4    0.920153
5    0.919519
Name: mae, dtype: float64

wqet_grader.grade("Project 3 Assessment", "Task 3.5.11", mae_series)

xxxxxxxxxx**Task 3.5.12:** Look through the results in `mae_series` and determine what value for `p` provides the best performance. Then build and train `final_model` using the best hyperparameter value.

best_model = AutoReg(y_train,lags=best_p).fit()

    "Project 3 Assessment", "Task 3.5.12", [isinstance(best_model.model, AutoReg)]

xxxxxxxxxx**Task 3.5.13:** Calculate the training residuals for `best_model` and assign the result to `y_train_resid`. **Note** that your name of your Series should be `"residuals"`.

y_train_resid.name = "residuals"

timestamp
2018-01-02 09:00:00+03:00   -0.560158
2018-01-02 10:00:00+03:00   -2.193105
2018-01-02 11:00:00+03:00   -0.026408
2018-01-02 12:00:00+03:00    0.820217
2018-01-02 13:00:00+03:00   -0.078009
Freq: H, Name: residuals, dtype: float64

wqet_grader.grade("Project 3 Assessment", "Task 3.5.13", y_train_resid.tail(1500))

xxxxxxxxxx**Task 3.5.14:** Create a histogram of `y_train_resid`. Be sure to label the x-axis as `"Residuals"` and the y-axis as `"Frequency"`. Use the title `"Best Model, Training Residuals"`.

plt.title("Best Model, Training Residuals")

    wqet_grader.grade("Project 3 Assessment", "Task 3.5.14", file)

xxxxxxxxxx**Task 3.5.15:** Create an ACF plot for `y_train_resid`. Be sure to label the x-axis as `"Lag [hours]"` and y-axis as `"Correlation Coefficient"`. Use the title `"Dar es Salaam, Training Residuals ACF"`.

plt.title("Dar es Salaam, Training Residuals ACF")

    wqet_grader.grade("Project 3 Assessment", "Task 3.5.15", file)

xxxxxxxxxx## Evaluate

xxxxxxxxxx**Task 3.5.16:** Perform walk-forward validation for your model for the entire test set `y_test`. Store your model's predictions in the Series `y_pred_wfv`. Make sure the name of your Series is `"prediction"` and the name of your Series index is `"timestamp"`.

    history = history.append(y_test[next_pred.index])

timestamp
2018-03-06 00:00:00+03:00    8.049707
2018-03-06 01:00:00+03:00    8.676565
2018-03-06 02:00:00+03:00    6.271593
2018-03-06 03:00:00+03:00    6.319877
2018-03-06 04:00:00+03:00    7.163557
Freq: H, Name: prediction, dtype: float64

wqet_grader.grade("Project 3 Assessment", "Task 3.5.16", y_pred_wfv)

xxxxxxxxxx**Task 3.5.17:** Submit your walk-forward validation predictions to the grader to see test mean absolute error for your model.

wqet_grader.grade("Project 3 Assessment", "Task 3.5.17", y_pred_wfv)

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In [40], line 1
----> 1 wqet_grader.grade("Project 3 Assessment", "Task 3.5.17", y_pred_wfv)

File /opt/conda/lib/python3.9/site-packages/wqet_grader/__init__.py:182, in grade(assessment_id, question_id, submission)
    177 def grade(assessment_id, question_id, submission):
    178   submission_object = {
    179     'type': 'simple',
    180     'argument': [submission]
    181   }
--> 182   return show_score(grade_submission(assessment_id, question_id, submission_object))

File /opt/conda/lib/python3.9/site-packages/wqet_grader/transport.py:146, in grade_submission(assessment_id, question_id, submission_object)
    144     raise Exception('Grader raised error: {}'.format(error['message']))
    145   else:
--> 146     raise Exception('Could not grade submission: {}'.format(error['message']))
    147 result = envelope['data']['result']
    149 # Used only in testing

Exception: Could not grade submission: Could not verify access to this assessment: Received error from WQET submission API: You have already passed this course!

xxxxxxxxxx# Communicate Results

xxxxxxxxxx**Task 3.5.18:** Put the values for `y_test` and `y_pred_wfv` into the DataFrame `df_pred_test` (don't forget the index). Then plot `df_pred_test` using plotly express. Be sure to label the x-axis as `"Date"` and the y-axis as `"PM2.5 Level"`. Use the title `"Dar es Salaam, WFV Predictions"`.

fig.write_image("images/3-5-18.png", scale=1, height=500, width=700)

    wqet_grader.grade("Project 3 Assessment", "Task 3.5.18", file)

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In [43], line 2
      1 with open("images/3-5-18.png", "rb") as file:
----> 2     wqet_grader.grade("Project 3 Assessment", "Task 3.5.18", file)

File /opt/conda/lib/python3.9/site-packages/wqet_grader/__init__.py:182, in grade(assessment_id, question_id, submission)
    177 def grade(assessment_id, question_id, submission):
    178   submission_object = {
    179     'type': 'simple',
    180     'argument': [submission]
    181   }
--> 182   return show_score(grade_submission(assessment_id, question_id, submission_object))

File /opt/conda/lib/python3.9/site-packages/wqet_grader/transport.py:146, in grade_submission(assessment_id, question_id, submission_object)
    144     raise Exception('Grader raised error: {}'.format(error['message']))
    145   else:
--> 146     raise Exception('Could not grade submission: {}'.format(error['message']))
    147 result = envelope['data']['result']
    149 # Used only in testing

Exception: Could not grade submission: Could not verify access to this assessment: Received error from WQET submission API: You have already passed this course!

xxxxxxxxxx---

​x

<font size="+3"><strong>Time Series: Core Concepts</strong></font>

from IPython.display import YouTubeVideo

# Model Types

## Autoregression Models

## ARMA Models

YouTubeVideo("voryLhxiPzE")

And in this video, Ritvik talks about the ARMA model we use in Project 3.

YouTubeVideo("HhvTlaN06AM")

# Plots

## ACF Plot

When we've worked with autocorrelations in the past, we've treated them like static relationships, but that's not always how they work. Sometimes, we want to actually see how those autocorrelations change over time, which means we need to think of them as *functions*. When we create a visual representation of an autocorrelation function (ACF), we're making an **ACF plot**.

## PACF Plot

# Statistical Concepts

## Walk-Forward Validation

## Parameters

## Hyperparameters

## Rolling Averages

---

---

xxxxxxxxxx​

xxxxxxxxxx<font size="+3"><strong>Databases: SQL</strong></font>

from IPython.display import YouTubeVideo

xxxxxxxxxx# Working with SQL Databases

xxxxxxxxxxA database is a collection of interrelated data. The primary goal of a database is to store and retrieve information in a convenient and efficient way. There are many types of databases. In this section, we will be dealing with a **relational database**. A relational database is a widely used database model that consists of a collection of uniquely named **tables** used to store information. The structure of a database model with its tables, constraints, and relationships is called a **schema**.

xxxxxxxxxx# Connecting to a Database

xxxxxxxxxx## ipython-sql

xxxxxxxxxx### Magic Commands

xxxxxxxxxxJupyter notebooks can run code that is not valid Python code but still affect the notebook . These special commands are called magic commands. Magic commands can have a range of properties. Some commonly used magic functions are below:

xxxxxxxxxx### ipython-sql

xxxxxxxxxx`ipython-sql` allows you to write SQL code directly in a Jupyter Notebook. The `%sql` (or `%%sql`) magic command is added to the beginning of a code block and then SQL code can be written.

xxxxxxxxxx### Connecting with ipython-sql

xxxxxxxxxxWe can connect to a database using the %sql magic function:

%sql sqlite:////home/jovyan/nepal.sqlite

xxxxxxxxxx## sqlite3

xxxxxxxxxxWe can also connect to the same database using the sqlite3 package:

conn = sqlite3.connect("/home/jovyan/nepal.sqlite")

xxxxxxxxxx# Querying a Database

xxxxxxxxxx## Building Blocks of the Basic Query

xxxxxxxxxxThere are six common clauses used for querying data:

xxxxxxxxxx## SELECT and FROM

xxxxxxxxxxYou can use `SELECT *` to select all columns in a table. `FROM` specifies the table in the database to query. `LIMIT 5` will select only the first five rows.

FROM id_map

xxxxxxxxxxYou can also use `SELECT` to select certain columns in a table

SELECT household_id,

xxxxxxxxxx<font size="+1">Practice</font>

%%sql

xxxxxxxxxxWe can also assign an **alias** or temporary name to a column using the `AS` command. Aliases can also be used on a table. See the example below, which assigns the alias `household_number` to `household_id`

SELECT household_id AS household_number,

xxxxxxxxxx<font size="+1">Practice</font>

%%sql

xxxxxxxxxx## Filtering and Sorting Data

xxxxxxxxxxSQL provides a variety of comparison operators that can be used with the WHERE clause to filter the data.

xxxxxxxxxxFor example, to select the first 5 homes in Ramechhap (district `2`):

%%sql

xxxxxxxxxx<font size="+1">Practice</font>

%%sql

xxxxxxxxxx## Aggregating Data

xxxxxxxxxxAggregation functions take a collection of values as inputs and return one value as the output. The table below gives the frequently used built-in aggregation functions:

xxxxxxxxxxUse the `COUNT` function to find the number of observations in the `id_map` table that come from Ramechhap (district `2`):

WHERE district_id = 2

xxxxxxxxxxAggregation functions are frequently used with a `GROUP BY` clause to perform the aggregation on groups of data. For example, the query below returns the count of observations in each District:

GROUP BY district_id

xxxxxxxxxx `DISTINCT` is a keyword to select unique records in a query result. For example, if we want to know the unique values in the `district_id` column:

SELECT distinct(district_id)

xxxxxxxxxx<font size="+1">Practice</font>

%%sql

xxxxxxxxxx`DISTINCT` and `COUNT` can be used in combination to count the number of distinct records. For example, if we want to know the number of unique values in the `district_id` column:

SELECT count(distinct(district_id))

xxxxxxxxxx<font size="+1">Practice</font>

%%sql

xxxxxxxxxx# Joining Tables

xxxxxxxxxxJoins link data from two or more tables together by using a column that is common between the two tables. The basic syntax for a join is below, where `table1` and `table2` refer to the two tables being joined, `column1` and `column2` refer to columns to be returned from both tables, and `ID` refers to the common column in the two tables.

xxxxxxxxxxWe'll explore the concept of joins by first identifying a single household that we'd like to pull in building information for. For example, let's say we want to see the corresponding `foundation_type` for the first home in Ramechhap (District 1). We'll start by looking at this single record in the `id_map` table.

WHERE district_id = 2

xxxxxxxxxxThis household has `building_id` equal to 23. Let's look at the `foundation_type` for this building, by filtering the `building_structure` table to find this building.

FROM building_structure

xxxxxxxxxxTo join the two tables and limit the results to `building_id = 23`:

JOIN building_structure ON id_map.building_id = building_structure.building_id

xxxxxxxxxxIn addition to the basic `JOIN` clause, specific join types can be specified, which specify whether the common column needs to be in one, both, or either of the two tables being joined. The different join types are below. The left table is the table specified first, immediately after the `FROM` clause and the right table is the table specified after the `JOIN` clause. If the generic `JOIN` clause is used, then by default the `INNER JOIN` will be used.

xxxxxxxxxxThe video below outlines the main types of joins:

YouTubeVideo("2HVMiPPuPIM")

xxxxxxxxxx<font size="+1">Practice</font>

%%sql

xxxxxxxxxx# Using pandas with SQL Databases

xxxxxxxxxxTo save the output of a query into a pandas DataFrame, we will use connect to the SQLite database using the SQLite3 package:

conn = sqlite3.connect("/home/jovyan/nepal.sqlite")

xxxxxxxxxxTo run a query using `sqlite3`, we need to store the query as a string. For example, the variable below called `query` is a string containing a query which returns the first 10 rows from the `id_map` table:

    FROM id_map

xxxxxxxxxxTo save the results of the query into a pandas DataFrame, use the `pd.read_sql()` function. The optional parameter `index_col` can be used to set the index to a specific column from the query.

df = pd.read_sql(query, conn, index_col="building_id")

xxxxxxxxxx<font size="+1">Practice</font>

query = ...

xxxxxxxxxx# References & Further Reading

xxxxxxxxxx---

xxxxxxxxxx​

xxxxxxxxxx<font size="+3"><strong>3.3. Autoregressive Models</strong></font>

xxxxxxxxxx import warnings​import matplotlib.pyplot as pltimport pandas as pdimport plotly.express as pxfrom IPython.display import VimeoVideofrom pymongo import MongoClientfrom sklearn.metrics import mean_absolute_errorfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacffrom statsmodels.tsa.ar_model import AutoReg​warnings.simplefilter(action="ignore", category=FutureWarning)

/opt/conda/lib/python3.9/site-packages/statsmodels/tsa/base/tsa_model.py:7: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  from pandas import (to_datetime, Int64Index, DatetimeIndex, Period,
/opt/conda/lib/python3.9/site-packages/statsmodels/tsa/base/tsa_model.py:7: FutureWarning: pandas.Float64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  from pandas import (to_datetime, Int64Index, DatetimeIndex, Period,

xxxxxxxxxx VimeoVideo("665851858", h="e39fc3d260", width=600)

xxxxxxxxxx# Prepare Data

xxxxxxxxxx## Import

xxxxxxxxxx VimeoVideo("665851852", h="16aa0a56e6", width=600)

xxxxxxxxxx**Task 3.3.1:** Complete to the create a client to connect to the MongoDB server, assigns the `"air-quality"` database to `db`, and assigned the `"nairobi"` connection to `nairobi`.

xxxxxxxxxx client = MongoClient(host="localhost", port=27017)db = client["air-quality"]​nairobi = db["nairobi"]

xxxxxxxxxx VimeoVideo("665851840", h="e048425f49", width=600)

xxxxxxxxxx**Task 3.3.2:** Change the `wrangle` function below so that it returns a Series of the resampled data instead of a DataFrame.

xxxxxxxxxx def wrangle(collection):    results = collection.find(        {"metadata.site": 29, "metadata.measurement": "P2"},        projection={"P2": 1, "timestamp": 1, "_id": 0},    )​    # Read data into DataFrame    df = pd.DataFrame(list(results)).set_index("timestamp")​    # Localize timezone    df.index = df.index.tz_localize("UTC").tz_convert("Africa/Nairobi")​    # Remove outliers    df = df[df["P2"] < 500]​    # Resample to 1hr window    y = df["P2"].resample("1H").mean().fillna(method='ffill')​    return y

xxxxxxxxxx**Task 3.3.3:** Use your wrangle function to read the data from the `nairobi` collection into the Series `y`.

xxxxxxxxxx y = wrangle(nairobi)y.head()

timestamp
2018-09-01 03:00:00+03:00    17.541667
2018-09-01 04:00:00+03:00    15.800000
2018-09-01 05:00:00+03:00    11.420000
2018-09-01 06:00:00+03:00    11.614167
2018-09-01 07:00:00+03:00    17.665000
Freq: H, Name: P2, dtype: float64

xxxxxxxxxx # Check your workassert isinstance(y, pd.Series), f"`y` should be a Series, not type {type(y)}"assert len(y) == 2928, f"`y` should have 2928 observations, not {len(y)}"assert y.isnull().sum() == 0

xxxxxxxxxx## Explore

xxxxxxxxxx VimeoVideo("665851830", h="85f58bc92b", width=600)

xxxxxxxxxx #y.corr(y) #correlation between y it self is 1#y.corr(y.shift(1)) #correlation between y its one hour lag value is 0.65#y.corr(y.shift(2)) #correlation between y its one hour lag value is 0.4#y.corr(y.shift(3)) #correlation between y its one hour lag value is 0.3#we can find the same correlations using ACF plot

xxxxxxxxxx**Task 3.3.4:** Create an ACF plot for the data in `y`. Be sure to label the x-axis as `"Lag [hours]"` and the y-axis as `"Correlation Coefficient"`.

xxxxxxxxxx fig, ax = plt.subplots(figsize=(15, 6))plot_acf(y,ax=ax)plt.xlabel("Lag [hours]")plt.ylabel("Correlation Coefficient");

xxxxxxxxxx VimeoVideo("665851811", h="ee3a2b5c24", width=600)

xxxxxxxxxx**Task 3.3.5:** Create an PACF plot for the data in `y`. Be sure to label the x-axis as `"Lag [hours]"` and the y-axis as `"Correlation Coefficient"`.

xxxxxxxxxx fig, ax = plt.subplots(figsize=(15, 6))plot_pacf(y,ax=ax)plt.xlabel("Lag [hours]")plt.ylabel("Correlation Coefficient");

xxxxxxxxxx## Split

xxxxxxxxxx VimeoVideo("665851798", h="6c191cd94c", width=600)

xxxxxxxxxx**Task 3.3.6:** Split `y` into training and test sets. The first 95% of the data should be in your training set. The remaining 5% should be in the test set.

xxxxxxxxxx cutoff_test = int(len(y)*0.95)​y_train = y.iloc[:cutoff_test]y_test = y.iloc[cutoff_test:]

xxxxxxxxxx len(y_train)+len(y_test)==len(y)

True

xxxxxxxxxx# Build Model

xxxxxxxxxx## Baseline

xxxxxxxxxx**Task 3.3.7:** Calculate the baseline mean absolute error for your model.

xxxxxxxxxx y_train_mean = y_train.mean()y_pred_baseline = [y_train_mean] * len(y_train)mae_baseline = mean_absolute_error(y_train, y_pred_baseline)​print("Mean P2 Reading:", round(y_train_mean, 2))print("Baseline MAE:", round(mae_baseline, 2))

Mean P2 Reading: 9.22
Baseline MAE: 3.71

xxxxxxxxxx## Iterate

xxxxxxxxxx VimeoVideo("665851769", h="94a4296cde", width=600)

xxxxxxxxxx**Task 3.3.8:** Instantiate an [`AutoReg`](https://www.statsmodels.org/stable/generated/statsmodels.tsa.ar_model.AutoReg.html) model and fit it to the training data `y_train`. Be sure to set the `lags` argument to `26`.

xxxxxxxxxx model = AutoReg(y_train,lags=26).fit()

xxxxxxxxxx VimeoVideo("665851746", h="1a4511e883", width=600)

xxxxxxxxxx**Task 3.3.9:** Generate a list of training predictions for your model and use them to calculate your training mean absolute error.

xxxxxxxxxx y_pred = model.predict().dropna()training_mae = mean_absolute_error(y_train.iloc[26:],y_pred)print("Training MAE:", training_mae)

Training MAE: 2.2809871656467036

xxxxxxxxxx VimeoVideo("665851744", h="60d053b455", width=600)

xxxxxxxxxx**Task 3.3.10:** Use `y_train` and `y_pred` to calculate the residuals for your model.

xxxxxxxxxx #y_train_resid = y_train-y_pred #which is equal to or same results with below codey_train_resid = model.residy_train_resid.tail()

timestamp
2018-12-25 19:00:00+03:00   -0.392002
2018-12-25 20:00:00+03:00   -1.573180
2018-12-25 21:00:00+03:00   -0.735747
2018-12-25 22:00:00+03:00   -2.022221
2018-12-25 23:00:00+03:00   -0.061916
Freq: H, dtype: float64

xxxxxxxxxx VimeoVideo("665851712", h="9ff0cdba9c", width=600)

xxxxxxxxxx**Task 3.3.11:** Create a plot of `y_train_resid`.

xxxxxxxxxx fig, ax = plt.subplots(figsize=(15, 6))y_train_resid.plot(ylabel="Residual Value",ax=ax)

<AxesSubplot:xlabel='timestamp', ylabel='Residual Value'>

xxxxxxxxxx VimeoVideo("665851702", h="b494adc297", width=600)

xxxxxxxxxx**Task 3.3.12:** Create a histogram of `y_train_resid`.

xxxxxxxxxx y_train_resid.hist()plt.xlabel("Residual Value")plt.ylabel("Frequency")plt.title("AR(26),Distribution of Residuals")

Text(0.5, 1.0, 'AR(26),Distribution of Residuals')

xxxxxxxxxx VimeoVideo("665851684", h="d6d782a1f3", width=600)

xxxxxxxxxx**Task 3.3.13:** Create an ACF plot of `y_train_resid`.

xxxxxxxxxx fig, ax = plt.subplots(figsize=(15, 6))plot_acf(y_train_resid,ax=ax)

xxxxxxxxxx## Evaluate

xxxxxxxxxx VimeoVideo("665851662", h="72e767e121", width=600)

xxxxxxxxxx # #method showed in video to work on test data and check mean_absolute_error# y_test.index.min()#to get the first value in test data using index# y_test.index.max()##to get the first value in test data using index# y_pred_test = model.predict(y_test.index.min(),y_test.index.max()).dropna()# test_mae = mean_absolute_error(y_test,y_pred_test)# print("Test MAE:", test_mae)

xxxxxxxxxx #method used in training modelmodel = AutoReg(y_test,lags=26).fit()

xxxxxxxxxx**Task 3.3.14:** Calculate the test mean absolute error for your model.

xxxxxxxxxx y_pred_test = model.predict().dropna()test_mae = mean_absolute_error(y_test.iloc[26:],y_pred_test)print("Test MAE:", test_mae)

Test MAE: 1.1581601912057071

xxxxxxxxxx**Task 3.3.15:** Create a DataFrame `test_predictions` that has two columns: `"y_test"` and `"y_pred"`. The first should contain the true values for your test set, and the second should contain your model's predictions. Be sure the index of `test_predictions` matches the index of `y_test`.

xxxxxxxxxx df_pred_test = pd.DataFrame(    {"y_test": y_test, "y_pred": y_pred_test}, index=y_test.index)

xxxxxxxxxx VimeoVideo("665851628", h="29b43e482e", width=600)

xxxxxxxxxx**Task 3.3.16:** Create a time series plot for the values in `test_predictions` using plotly express. Be sure that the y-axis is properly labeled as `"P2"`.

xxxxxxxxxx fig = px.line(df_pred_test, labels={"value": "P2"})fig.show()

xxxxxxxxxx VimeoVideo("665851599", h="bb30d96e43", width=600)

xxxxxxxxxx**Task 3.3.17:** Perform walk-forward validation for your model for the entire test set `y_test`. Store your model's predictions in the Series `y_pred_wfv`.

xxxxxxxxxx %%capture​y_pred_wfv = pd.Series()history = y_train.copy()for i in range(len(y_test)):    model = AutoReg(history,lags=26).fit()    next_pred = model.forecast()    y_pred_wfv = y_pred_wfv.append(next_pred)    history = history.append(y_test[next_pred.index])    ​

xxxxxxxxxx len(y_pred_wfv)==len(y_test)

True

xxxxxxxxxx history.tail(1)

2019-01-01 02:00:00+03:00    18.803333
Freq: H, Name: P2, dtype: float64

xxxxxxxxxx ​

xxxxxxxxxx model = AutoReg(history,lags=26).fit()

xxxxxxxxxx model.forecast()

2019-01-01 03:00:00+03:00    14.110356
Freq: H, dtype: float64

xxxxxxxxxx y_test.head(1)

timestamp
2018-12-26 00:00:00+03:00    5.679091
Freq: H, Name: P2, dtype: float64

xxxxxxxxxx ​

xxxxxxxxxx VimeoVideo("665851568", h="a764ab5416", width=600)

xxxxxxxxxx**Task 3.3.18:** Calculate the test mean absolute error for your model.

xxxxxxxxxx test_mae = mean_absolute_error(y_test,y_pred_wfv)print("Test MAE (walk forward validation):", round(test_mae, 2))

Test MAE (walk forward validation): 1.4

xxxxxxxxxx# Communicate Results

xxxxxxxxxx VimeoVideo("665851553", h="46338036cc", width=600)

xxxxxxxxxx**Task 3.3.19:** Print out the parameters for your trained model.

xxxxxxxxxx #unlike in scikit learn finding intercept and coefficients in statsmodels is very easy. below code will give those intercept and coefficient valuesprint(model.params)

intercept    2.021763
P2.L1        0.587991
P2.L2        0.019749
P2.L3        0.023426
P2.L4        0.026807
P2.L5        0.044019
P2.L6       -0.101889
P2.L7        0.029814
P2.L8        0.049894
P2.L9       -0.017513
P2.L10       0.032585
P2.L11       0.064034
P2.L12       0.005978
P2.L13       0.018366
P2.L14      -0.007634
P2.L15      -0.016571
P2.L16      -0.015587
P2.L17      -0.035295
P2.L18       0.000929
P2.L19      -0.003956
P2.L20      -0.020668
P2.L21      -0.012730
P2.L22       0.052431
P2.L23       0.074405
P2.L24      -0.024263
P2.L25       0.090515
P2.L26      -0.088646
dtype: float64

xxxxxxxxxx VimeoVideo("665851529", h="39284d37ac", width=600)

xxxxxxxxxx**Task 3.3.20:** Put the values for `y_test` and `y_pred_wfv` into the DataFrame `df_pred_test` (don't forget the index). Then plot `df_pred_test` using plotly express.

xxxxxxxxxx df_pred_test = pd.DataFrame(    {"y_test":y_test,"y_pred_wfv":y_pred_wfv}​)fig = px.line(df_pred_test,labels={"Value":"PM2.5"})fig.show()

xxxxxxxxxx---

xxxxxxxxxx​

xxxxxxxxxx<font size="+3"><strong>3.2. Linear Regression with Time Series Data</strong></font>

from sklearn.linear_model import LinearRegression

VimeoVideo("665412117", h="c39a50bd58", width=600)

xxxxxxxxxx# Prepare Data

xxxxxxxxxx## Import

VimeoVideo("665412469", h="135f32c7da", width=600)

xxxxxxxxxx**Task 3.2.1:** Complete to the create a client to connect to the MongoDB server, assign the `"air-quality"` database to `db`, and assign the `"nairobi"` connection to `nairobi`.

client = MongoClient(host="localhost",port=27017)

VimeoVideo("665412480", h="c20ed3e570", width=600)

xxxxxxxxxx**Task 3.2.2:** Complete the `wrangle` function below so that the `results` from the database query are read into the DataFrame `df`. Be sure that the index of `df` is the `"timestamp"` from the results.

    df.index = df.index.tz_localize("UTC").tz_convert("Africa/Nairobi")

VimeoVideo("665412496", h="d757475f7c", width=600)

xxxxxxxxxx**Task 3.2.3:** Use your wrangle function to read the data from the `nairobi` collection into the DataFrame `df`.

df = wrangle(nairobi)

assert any([isinstance(df, pd.DataFrame), isinstance(df, pd.Series)])

VimeoVideo("665412520", h="e03eefff07", width=600)

xxxxxxxxxx**Task 3.2.4:** Add to your `wrangle` function so that the `DatetimeIndex` for `df` is localized to the correct timezone, `"Africa/Nairobi"`. Don't forget to re-run all the cells above after you change the function.

# df.index.tz_localize("UTC").tz_convert("Africa/Nairobi")

assert df.index.tzinfo == pytz.timezone("Africa/Nairobi")

xxxxxxxxxx## Explore

VimeoVideo("665412546", h="97792cb982", width=600)

xxxxxxxxxx**Task 3.2.5:** Create a boxplot of the `"P2"` readings in `df`.

df["P2"].plot(kind="box",vert=False,ax=ax)

<AxesSubplot:>

VimeoVideo("665412573", h="b46049021b", width=600)

xxxxxxxxxx**Task 3.2.6:** Add to your `wrangle` function so that all `"P2"` readings above 500 are dropped from the dataset. Don't forget to re-run all the cells above after you change the function.

assert len(df) <= 32906

VimeoVideo("665412594", h="e56c2f6839", width=600)

xxxxxxxxxx**Task 3.2.7:** Create a time series plot of the `"P2"` readings in `df`.

fig,ax = plt.subplots(figsize=(15, 6))

VimeoVideo("665412601", h="a16c5a73fc", width=600)

xxxxxxxxxx**Task 3.2.8:** Add to your `wrangle` function to resample `df` to provide the mean `"P2"` reading for each hour. Use a forward fill to impute any missing values. Don't forget to re-run all the cells above after you change the function.

# df["P2"].resample("1H").mean().fillna(method="ffill").to_frame() #.isnull().sum()

assert len(df) <= 2928

VimeoVideo("665412649", h="d2e99d2e75", width=600)

xxxxxxxxxx**Task 3.2.9:** Plot the rolling average of the `"P2"` readings in `df`. Use a window size of `168` (the number of hours in a week).

df["P2"].rolling(168).mean().plot(ax=ax);

VimeoVideo("665412693", h="c3bca16aff", width=600)

xxxxxxxxxx**Task 3.2.10:** Add to your `wrangle` function to create a column called `"P2.L1"` that contains the mean`"P2"` reading from the previous hour. Since this new feature will create `NaN` values in your DataFrame, be sure to also drop null rows from `df`.

assert len(df) <= 11686

VimeoVideo("665412732", h="059e4088c5", width=600)

xxxxxxxxxx**Task 3.2.11:** Create a correlation matrix for `df`.

df.corr()

VimeoVideo("665412741", h="7439cb107c", width=600)

xxxxxxxxxx**Task 3.2.12:** Create a scatter plot that shows PM 2.5 mean reading for each our as a function of the mean reading from the previous hour. In other words, `"P2.L1"` should be on the x-axis, and `"P2"` should be on the y-axis. Don't forget to label your axes!

ax.plot([0,120],[0,120],linestyle="--",color="orange")

[<matplotlib.lines.Line2D at 0x7f7eae8ce940>]

xxxxxxxxxx## Split

VimeoVideo("665412762", h="a5eba496f7", width=600)

xxxxxxxxxx**Task 3.2.13:** Split the DataFrame `df` into the feature matrix `X` and the target vector `y`. Your target is `"P2"`.

feature=["P2.L1"]

VimeoVideo("665412785", h="03118eda71", width=600)

xxxxxxxxxx**Task 3.2.14:** Split `X` and `y` into training and test sets. The first 80% of the data should be in your training set. The remaining 20% should be in the test set.

X_train, y_train = X.iloc[:cutoff],y.iloc[:cutoff]

len(X_train)+len(X_test)==len(X)

True

xxxxxxxxxx# Build Model

xxxxxxxxxx## Baseline

xxxxxxxxxx**Task 3.2.15:** Calculate the baseline mean absolute error for your model.

mae_baseline = mean_absolute_error(y_train,y_pred_baseline)

Mean P2 Reading: 9.27
Baseline MAE: 3.89

xxxxxxxxxx## Iterate

xxxxxxxxxx**Task 3.2.16:** Instantiate a [`LinearRegression`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html) model named `model`, and fit it to your training data.

model = LinearRegression()

LinearRegression()

LinearRegression()

xxxxxxxxxx## Evaluate

VimeoVideo("665412844", h="129865775d", width=600)

xxxxxxxxxx**Task 3.2.17:** Calculate the training and test mean absolute error for your model.

training_mae = mean_absolute_error(y_train,model.predict(X_train))

Training MAE: 2.46
Test MAE: 1.8

xxxxxxxxxx# Communicate Results

xxxxxxxxxx**Task 3.2.18:** Extract the intercept and coefficient from your `model`.

print(f"P2 = {intercept} + ({coefficient} * P2.L1)")

P2 = 3.36 + (0.64 * P2.L1)

VimeoVideo("665412870", h="318d69683e", width=600)

xxxxxxxxxx**Task 3.2.19:** Create a DataFrame `df_pred_test` that has two columns: `"y_test"` and `"y_pred"`. The first should contain the true values for your test set, and the second should contain your model's predictions. Be sure the index of `df_pred_test` matches the index of `y_test`.

df_pred_test = pd.DataFrame({"y_test":y_test,"y_pred":model.predict(X_test)})

VimeoVideo("665412891", h="39d7356a26", width=600)

xxxxxxxxxx**Task 3.2.20:** Create a time series line plot for the values in `test_predictions` using plotly express. Be sure that the y-axis is properly labeled as `"P2"`.

fig = px.line(df_pred_test)

xxxxxxxxxx---

Rolling Months	Rolling sum sales
2022-01,02,03	55
2022-02,03,04	60
2022-03,04,05	60
2022-04,05,06	65

	P2
timestamp
2018-09-01 00:00:02.472	34.43
2018-09-01 00:05:03.941	30.53
2018-09-01 00:10:04.374	22.80
2018-09-01 00:15:04.245	13.30
2018-09-01 00:20:04.869	16.57

Magic Command	Description of Command
`%pwd`	Print the current working directory
`%cd`	Change the current working directory
`%ls`	List the contents of the current directory
`%history`	Show the history of the In [ ]: commands

Clause Name	Definition
`SELECT`	Determines which columns to include in the query's result
`FROM`	Identifies the table from which to query the data from
`WHERE`	filters data
`GROUP BY`	groups rows by common values in columns
`HAVING`	filters out unwanted groups from GROUP BY
`ORDER BY`	Orders the rows using one or more columns
`LIMIT`	Outputs the specified number of rows

Comparison Operator	Description
=	Equal
>	Greater than
<	Less than
>=	Greater than or equal to
<=	Less than or equal to
<> or !=	Not equal to
LIKE	String comparison test

	P2
timestamp
2018-10-01 00:04:20.554	0.01
2018-10-01 00:07:47.504	0.01
2018-10-01 00:11:14.382	0.01
2018-10-01 00:14:41.239	0.01
2018-10-01 00:18:05.938	0.01

	0	8	16	24
0	4.1715	2.9384	2.9201	2.9144
1	3.3506	2.9149	2.9294	2.9136
2	3.1057	2.9080	2.9136	2.8979

Aggregation Function	Definition
`MIN`	Return the minimum value
`MAX`	Return the largest value
`SUM`	Return the sum of values
`AVG`	Return the average of values
`COUNT`	Return the number of observations

JOIN Type	Definition
`INNER JOIN`	Returns rows where ID is in both tables
`LEFT JOIN`	Returns rows where ID is in the left table. Return NA for values in column, if ID is not in right table.
`RIGHT JOIN`	Returns rows where ID is in the right table. Return NA for values in column, if ID is not in left table.
`FULL JOIN`	Returns rows where ID is in either table. Return NA for values in column, if ID is not in either table.
WQU WorldQuant University Applied Data Science Lab QQQQ

	P2	P2.L1
timestamp
2018-09-01 04:00:00+03:00	15.800000	17.541667
2018-09-01 05:00:00+03:00	11.420000	15.800000
2018-09-01 06:00:00+03:00	11.614167	11.420000
2018-09-01 07:00:00+03:00	17.665000	11.614167
2018-09-01 08:00:00+03:00	21.016667	17.665000
2018-09-01 09:00:00+03:00	22.589167	21.016667
2018-09-01 10:00:00+03:00	18.605833	22.589167
2018-09-01 11:00:00+03:00	14.022500	18.605833
2018-09-01 12:00:00+03:00	13.150000	14.022500
2018-09-01 13:00:00+03:00	12.806667	13.150000

	y_test	y_pred
timestamp
2018-12-07 17:00:00+03:00	7.070000	8.478927
2018-12-07 18:00:00+03:00	8.968333	7.865485
2018-12-07 19:00:00+03:00	11.630833	9.076421
2018-12-07 20:00:00+03:00	11.525833	10.774814
2018-12-07 21:00:00+03:00	9.533333	10.707836

Month	sales
2022-01	10
2022-02	20
2022-03	25
2022-04	15
2022-05	20
2022-06	30

	P2	P2.L1
P2	1.000000	0.650679
P2.L1	0.650679	1.000000

Open Tabs

Kernels

Terminals

Autoregression¶

Cleaning the Data¶

Localizing the Timezone¶

Resampling Data¶

Adding a Lag¶

Exploring the Data¶

Time Series Line Plot¶

Creating an ACF Plot¶

Creating a PACF Plot¶

Working with Rolling Windows¶

Splitting the Data in pandas¶

Building the Model¶

Baseline¶

Iterating¶

Residuals¶

ACF Plots¶

Evaluating the Model¶

Out-of-Sample Predictions¶

Walk-forward Validation¶

Communicating the Results¶

ARMA Models & Hyperparameters¶

Cleaning the Data¶

Exploring the Data¶

Splitting the Data¶

Building the Model¶

Baseline¶

Iterating¶

Hyperparameters¶

Communicating the Results¶

ARCH and GARCH Models¶

Calculating Returns¶

Building the Model¶

Common Metrics¶

Standardized Residuals¶

Evaluation¶

Forecasting¶

References & Further Reading¶

Prepare Data¶

Connect¶

Explore¶

Import¶

Working with PyMongo¶

Databases¶

Servers and Clients¶

Semi-structured Data¶

Exploring a Database¶

Collections¶

Documents¶

Retrieving Data¶

Key-Value Pairs¶

Metadata¶

Identifiers¶

Analyzing Data¶

Updating Documents¶

Aggregation¶

Finishing Up¶

References & Further Reading¶

Prepare Data¶

Import¶

Explore¶

Split¶

Build Model¶

Baseline¶

Iterate¶

Evaluate¶

Communicate Results¶

Calculate Summary Statistics for a DataFrame or Series¶

Select a Series from a DataFrame¶

Subset a DataFrame by Selecting One or More Columns¶

Subset the Columns of a DataFrame Based on Data Types¶

Working with value_counts in a Series¶

Series and Groupby¶

Subset a DataFrame's Columns Based on the Column Data Types¶

Subset a DataFrame's columns based on column names¶

Pivot Tables¶

Subsetting with Masks¶

Reshape a DataFrame based on column values¶

Working with `value_counts` in a Series¶

Series and `Groupby`¶